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depths up to 18 km w.e. in different rocks and in water are calculated. The results are particularly 
collated with a great body of the ground-level, underground, and underwater muon data. In the 
hadron-cascade calculations, we take into account the logarithmic growth with energy of inelastic 
cross sections and pion, kaon, and nucleon generation in pion-nucleus collisions. For evaluating the 
prompt muon contribution to the muon flux, we apply the two phenomenological approaches to 
the charm production problem: the recombination quark-parton model and the quark-gluon string 
model. We give simple fitting formulas describing our numerical results. To solve the muon transport 
equation at large depths of a homogeneous medium, we use a semianalytical method, which allows 
the inclusion of an arbitrary (decreasing) muon spectrum at the medium boundary and real energy 
dependence of muon energy losses. Our analysis shows that at the depths up to 6-7 km w.e., 
essentially all underground data on the muon flux correlate with each other and with the predicted 
one for conventional (tt, K) muons, to within 10 %. However, the high-energy sea-level muon data 
as well as the data at high depths are contradictory and cannot be quantitatively described by a 
single nuclear- cascade model. 
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I. INTRODUCTION 



The flux of cosmic-ray muons in the atmosphere, underground, and underwater provides a way of testing the 
inputs of nuclear cascade models, that is, parameters of the primary cosmic-ray flux (energy spectrum, chemical 
composition) and particle interactions at high energies. In particular, measurements of the muon energy spectra, 
angular distributions and the depth-intensity relation (DIR) have much potential for yielding information about the 
mechanism of charm production in hadron-nucleus collisions at energies beyond the reach of accelerator experiments. 
This information is a subject of great current interest for particle physics and yet is a prime necessity in high- 
energy and very high-energy neutrino astronomy |^ . Indeed, the basic and unavoidable background for many future 
astrophysical experiments with full-size underwater or underice neutrino telescopes will be an effect of the atmospheric 
neutrino flux of energies from about 1 TeV to tens of PcV. However, in the absence of a generally recognized and 
tried model for charm hadroproduction, the current estimates of the z^^ and (most notably) Ve backgrounds have 
inadmissibly wide scatter even at multi-TcV neutrino energies, which shoots up with energy. At E^j ~ 100 TeV, 
different estimates of the and Ve spectra vary within a few orders of magnitude (see Refs. for reviews and 

references) . 

The present state of the art of predicting the atmospheric neutrino flux seems to be more satisfactory at energies 
below a few TeV. However, the theory meets more rigid requirements on accuracy of the calculations here Q: for 
an unambiguous treatment of the current data on up-going (atmospheric neutrino induced) muon flux, the neutrino 
flux must be calculated with a 10 % accuracy at least, whereas the uncertainties in the required input data (primary 
spectrum, cross sections for light meson production, etc.) hinder to gain these ends. Because of this, a vital question 
is a normalization of the calculated (model-dependent) atmospheric neutrino flux and the muon flux is perhaps the 
only tool for such a normalization. The point is that atmospheric muons and neutrinos are generated in just the same 
processes. Therefore, the accuracy of the neutrino flux calculation can be improved by forcing the poorly known input 
parameters of the cascade model to flt the data on the muon flux. 

The sea-level muon data obtained by direct measurements with magnetic spectrometers are crucial but still insuf- 
flcient for this purpose. The fact is that numerous sea-level measurements (see e.g. Refs. ||7|"|l7| for the vertical muon 
flux, Refs. p^Jl^ for near-horizontal flux, and Ref. |20i for a compilation of the data) are in rather poor agreement to 
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one another, even though each of the experiments by itself typically has very good statistical accuracy. This is true 
to a greater or lesser extent everywhere over the whole energy region accessible to the ground-based installations. 

On the other hand, a quite representative array of data on cosmic-ray muon DIR in rock and, to a lesser extent, 
in water has been accumulated. Underground muon experiments may number in the tens in a span of sixty years 
(see Refs. and also for reviews and further references). It should be noted that the results of many 

early measurements, specifically those performed at shallow depths, have not lost their significance today, considering 
that modern experiments principally aim at greater depths. Underwater muon experiments have over 30 years of 
history [pT|-[59|] and it is believed that they will gain in importance with the progress of high-energy neutrino telescopes. 

It may be somewhat unexpected but the underground data are more self-consistent in comparison with ground-level 
data, at least for depths to about 6 km w.e. (corresponding roughly to 3-4 TeV of muon energy at sea level) and 
hence they provide a useful check on nuclear cascade models. There is a need to piece together all these data in order 
to extract some physical inferences thereof. Also, it would be useful to correlate the underground and underwater 
data with the results of the mentioned direct measurements of the sea-level muon spectrum |7|-p7| as well as with the 
data deduced by indirect routes ^3 36 3^Q,^ -66|. 

It is the purpose of this paper to discuss the above-mentioned data on the vertical muon flux at sea level, under- 
ground, and underwater in the context of a single calculation, with emphasis on the prompt muon problem. The 
implementation of the results to the normalization of the high-energy atmospheric neutrino flux will be discussed 
elsewhere pq. 

In Section y we discuss the employed model for the primary spectrum and composition as well as the nuclear- 
cascade model for production and propagation of high-energy nucleons, pions, and kaons in the atmosphere. Some 
required formulas for the atmospheric muon flux are given in Section III ; at the end of that Section, we give a simple 
parametrization for the calculated vertical spectrum of conventional (tt, K) muons at sea level. The models for charm 
hadroproduction, those are used in the present work to make an estimate of the prompt-muon (PM) contribution, 
are the concern of Section IV; the recombination quark-parton model is considered with some details. At the end of 
this Section, we present simple parametrizations for the predicted differential and integral PM spectra. In Section ^ 
we compare our predi ctio ns for the vertical muon spectra (differential and integral) with the direct and indirect data 
at sea level. Section 
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is concerned with muon propagation through matter. Calculation of the muon intensity 
at large depths is a rather nontrivial problem even though the muon energy spectrum at the medium boundary is 
assumed to be known; we briefly sketch our semianalytical approach to that problem. The comparison between the 
calculated muon DIR and the aforecited underground and underwater data is fully considered in Section VII. In 
Appendix ^ we give the model formulas for the spectra of muons from inclusive semileptonic decays of a I? meson 
and Ac hyperon in the lab. frame. In Appendix ^ we present a summary for the differential cross sections of the 
muon-matter interactions (direct e^e" pair production, bremsstrahlung, photonuclear interaction) as well as (fo r 
completeness sake) Sternheimer's formula for ionization energy loss. Our conclusions are presented in Section VIII. 



II. NUCLEAR-CASCADE MODEL 



A. Primary spectrum and composition 

For energies above 1 TeV we use the semiempirical model for the integral primary spectrum proposed by Nikol'sky 
et al. Ipq] (from here on we will call it "NSU model" ) : 



F{> Eo) = FoEp ^ A + 5a^ 



(2.1) 



Here Eq is the energy per particle in GeV, Fq = 1.16 cm~^s~^sr~^, 7 = 1.62 (±0.03), and ae = 0.4. The 5as 
specify the region of the "knee" in the primary spectrum. We adopt 5p — Q x 10^^ and 5a>4: = 10^^. These values 
correspond to the hypothesis which attributes the change in the energy spectrum of the primaries at Eq > 10'^ TeV 
to photodesintegration of nuclei with pion photoproduction by photons with energy ~ 70 eV inside the cosmic ray 
sources. The chemical composition is given with the following values for Ba'- Bi — 0.40 (±0.03), B4 = 0.21 (±0.03), 
Bi5 = 0.14 (±0.03), B26 = 0.13 (±0.03), and B51 = 0.12 (±0.04) for the five standard groups of nuclei. The numerical 
values of A indicate the average atomic weights in the groups. The corresponding differential spectrum is given by 



dF 
dE^ 



^FoE^^-^'^TBa 
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Eo 
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xSaEo/A 



7(1 + SaEo/A) 



(2.2) 
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The NSU approximation has been deduced from an analysis of fluctuations in the relative number of electrons and 
muons in extensive air showers and corresponds to the data on absolute intensities of primary protons and various 
nuclei at energies Eq > 1, 10'^, 10^ TeV/particle, and also to the data on the shape of the integral spectrum in the 
vicinity of the knee (see Ref . for specific sources of the data) . 

The model, on the whole, fits the modern data on the primary spectrum and composition from about 
100 GeV/particle up to 100 EeV/particle. Specifically, at Eq < 10^ TeV/particle, the model fits reasonably well 
the recent results of the COSMOS satellite experiment the JACEE balloon experiment and the BASJE air- 
shower experiment f7l|] . On the other hand, there is a strong discrepancy between the NSU model and the recent data 
of the Japan balloon-borne emulsion chamber experiment ||7^ , which indicates a milder knee shape than that found 
in the previous experiments, although the data of Ref. |7|] for the nuclear composition agree with the NSU model at 
Eq > 10 TeV/particle. The data for the spectrum and composition are most inconsistent in the vicinity of the knee 
[(10^ -f- 10^) TeV/particle] . Scanty experimental data favor a pure proton composition at Eq > 10^ TeV/particle rather 
than almost fixed composition predicted by the NSU model. In the connection it should be noted that an essential 
contribution to the deep underground flux of muons, in particular, ones originated from the decay of charmed hadrons 
(at depths below ~ 10 km w.e.), is given by the primaries with energies from the knee region. Thus the long-standing 
problem of the knee is closely allied to the PM problem. At the same time, the total intensity of underground/water 
muons is scarcely affected by the region Eq 3> 10* TeV/particle. Thus we will not discuss here the problem of the 
primary spectrum and composition at super- high energies (see Refs. |^,^ for current reviews). 



B. Nuclear cascade at high energies: Basic assumptions 



Our nuclear-cascade calculations at high energies are based on the analytical model of Ref. |Q which describes 
well all available experimental data on hadron spectra for various atmospheric depths and for energies from about 
1 TeV up to about 600 TeV. The processes of regeneration and overcharging of nucleons, and charged pions, as well 
as production of kaons, nucleons, and charmed particles in pion-nucleus collisions have been properly accounted for. 
Let us outline the basic assumptions of the model. 

(i) The nuclear component of the primary spectrum is replaced with a superposition of free nucleons. Eq. (2.2) 
transforming to the equivalent nucleon spectrum yields the following differential energy spectra of protons and neu- 
trons: 

^ = D°(B») = D,(E„) + IY1 

A>4 



A>4 



Here E^ is the nucleon energy (in GeV), 



■DaIEn) 



(1 + SaEn)'' 



sbSaE, 



N 



-f{l + SAEN) 



Vq = jBiFq = 0.75cm-2s-isr-i(GeV/nucleon)-i, and Ca = A^-'^Ba/Bi (A = 1, 4, 15, 26, 51). Outside the knee 
region we use the asymptotic formulas: 



Va{En) = 



CaVqE^^'<+'^ for£;w«£;«, 

(2) 
N ' 



1.25(5T"=CaI?oSv^^+"^'^ for En » E 



(2.3) 



where E 



(1) _ 

N 



6.5/Sa GeV/nucleon and E^^ 



0.6/6a GeV/nucleon. A numerical procedure is applied to smooth out 
the calculated spectra of secondary hadrons at energies around the knee region. 

(ii) We assume a logarithmic growth with energy of the total inelastic cross sections ctJ^°' for interactions of a hadron 
i with a nuclear target A. Such a dependence arises from a model for elastic amplitude of hadron-hadron collisions, 
based on the conception of double pomeron with the supercritical intercept | |7^ . For simplicity sake we will use also 
another consequence of this model: the asymptotic equality of the inelastic cross sections for any hadron. Thereby 
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(z = iV,^,i^,...) 



(2.4) 
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a,t E > El = 1 TeV. The following values of the parameters are adopted: aA = 19 mb, cr^^ = 275 mb {N = p,n), 
crO±^ = 212 mb, = 183 mb {K = K^,K°,K"). 

(iii) It is assumed that Feynman scaling holds in the fragmentation region of the inclusive processes iA fX, 
where i = p,n,TT^ , f — p,n,Tr^ , K^, K^, K^, and A is the "air nucleus" . So the normalized invariant inclusive cross 
sections (i?/cr^^) d^aiA^jx/d^P are energy independent at large x (where a; is the ratio of the final particle energy 
to that of the initial one). Let us denote 

Then the fractional moments ( "Z- factors" ) defined by 

Zf^{l) = / X^-^Wf^{x)dx (2.5) 

Jo 

are constant inside the regions with constant exponent 7 (that is outside the knee energy region in the primary 
spectrum). Table || shows fractional moments Zfi{'y) for the two values of 7 in the case where the incident particle i is 
a proton or 7r+ meson and f = p,n, n^, , K^. The moments for i = n and 7r~ can be derived using the well-known 
isotopic relations for the inclusive cross sections. To calculate the Z-factors for all reactions except iiA NX and 
ttA KX, we used a parametrization of ISR data put forward by Minorikawa and Mitsui ||7^. The quantities Z^-^ 
and Zkt^ were calculated from the two central moments, (x) and (x^), for the inclusive distributions obtained by 
Anisovich et al. [ f76[ in the framework of quasinuclear quark model. 

TABLE I. Fractional moments Zfi{"f) of inclusive distributions of nucleons, pious, and kaons for the two values of 7. 
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7 = 1.62 
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0.1990 


0.0763 


0.0474 


0.0318 


0.0067 


0.0023 


0.0045 




0.0070 


0.0060 


0.1500 


0.0552 


0.0120 


0.0120 


0.0100 










7 = 2.02 








P 


0.1980 


0.0585 


0.0257 


0.0162 


0.0039 


0.0012 


0.0026 


n+ 


0.0060 


0.0040 


0.1480 


0.0346 


0.0100 


0.0100 


0.0080 



(iv) The kaon regeneration (i.e. the processes K^A —> K^X, K^A K^X , etc.) is disregarded in our calculations. 
Also, we neglect the nucleon and pion production in kaon-nucleus collisions as well as pion production in kaon decays, 
which makes it possible to split up the total system of the transport equations into nucleon-pion part and kaon one. 
Our estimations show that the inclusion of the aforementioned effects will cause the muon flux to increase, but no 
more than by a few per cent. It is clear that similar effects for charmed particles are completely negligible. 

(v) At the stage of nuclear cascade (but, of course, not at the muon production stage) the decay of tt^ mesons (critical 
energy E" ~ 0.12 TeV) is neglected for directions close to vertical at pion energies > 1 TeV. This approximation 
greatly simplifies the description of the pion regeneration and the production of nucleons, kaons, and charmed particles 
in pion-nucleus collisions. 

C. Nucleon-pion cascade equations 

In line with the above-listed assumptions, the 4x4 system of transport equations for the nucleon-pion part of the 
cascade can be written 

'd_ 1 

dh ^ X,{E) 

(«, j — p, n, 7r+, TT^) with the boundary conditions 

Vj,(E,Q)=Vl{E), V^{E,0)^V"JE), (i?, 0) = P,- (i?, 0) = 0. 



V,{E,h) 



(2.6) 
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Here 'Di{E, h) is the differential energy spectrum of particles i at the atmospheric depth h, 

1 .n 1 



ME) 



A 



and -/Vo is the number of target nuc lei in 1 g of air. 

The solution to the system (2.6) can be found as an expansion in powers of the dimcnsionless parameter /iMa, 
where Xa — 1/(A^oO'a) — 14.5A5V- Within the power-behaved regions of the primary spectrum described by Eq. (^.3|), 
the solution is of the form 

'DpiE, /i) = i [N+{E, h) + N-{E, h)] , VniE, h)^^ [N+{E, h) - N-{E, h)] , 

iE,h)^^ [n+ {E,h) + n {E, h)] , v,-iE,h)^^[n+{E,h)-n-{E,h)], 



where 



V°(E) + K,V°(E)^ 
N^iE, h) = ' (^"^ + ^ ) 

Vl(E) + nVl (E) / A, \ ^ 



^W^iE) 



l + O 



Xa 



2j- 
1 



AS,-,(B) 



l + K'j'^iE) l-K'j'^{E) 



K-£{E) 2A%iE) 



2A^{E) 



(k, k = ±), 



riE) = Ji 



yK yK A 2 
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a; 



1 - z 



NN 



2X% 



2AS ■ 



■Jpp ■ 



kZ, 



np, 



1 ^ ^ z^ 

k1{E) X,{E) X\ ' 



7r+7r + 



The functions K'^'^{E) can be treated as the generalized absorption ranges. Not counting the processes of nucleon- 
antinucleon pair production by pions, the formulas for A'^^{E) are very simple: 



^%\{E) - A%{E), A%-{E) = A:{E), 



and thus 



N'-'{E,h) cx exp 



^%iE) 



IT'(i?, h) cx exp — -r— 7 — r — exp 



A-iE) 



The 0{h/XA) corrections were calculated in Ref. and it was demonstrated that they became important at 
h > 500 — 600 g/cm^. However these corrections are of no significance for present purposes, because the greater part 
of the atmospheric muon flux is generated on the depths h < 300 g/cm^. 



D. Kaon production and transport 



Kaon decay cannot be neglected even at very high energies; as a result the d iffere ntial energy spectra of kaons, 
VjiiE, h, "(9), depend on zenith angle i?. In line with approximation (iv) of Section II B and assuming isothermality of 
the atmosphere, the transport equation for kaons may be written as 



d_ 
dh 



1 



Aa-(S) 



gg(^) 
Eh 



VK{E,h,d)^GKiE,h), [K^K^.Kl] 



(2.7) 
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where E'j^{'d) = mKHosecd/TK is the kaon critical energy (at d < 75°), uik and tk are the kaon mass and hfetime, 
and Hq ~ 6.44 km is the parameter of the isothermal atmosphere. The source function GKiE,h) describes kaon 
production in NA and ttA c ollis ions. Taking into account the explicit form of the nucleon and pion spectra outside 
the knee region (see Section |ll C| ) , we have 



GK{E,h) 



\ JQ 



(2.8) 



where 



and 7^ = 7 + H/Xa- 

Upon integrating Eq. (2.7) with the source function ( |2.8| ) and neglecting the weak /i-dependence of the kaon Z- 
factors we obtain 



VK{E,h,'d) = / exp 



\k{E) 



r(£^f(i?))exp 



\k{E) 



E 



GK{E,h')dh' 

Zkn{i) (3^) N-K{E,h,§) + Z-K^{^) (^A J nj,(i?,/j,^) 



(2.9) 



nK(i?,/^,^) = ^.^(7) (^A^j E(-« h (^^^W' 



Here e/f (z?) = E'i^{d)/E + 1, F is the gamma- function, 7 is the incomplete gamma- function. 



Tie) 



and 



K-£-{E) A-^;{E) \k{E) 

The approximate solution (|2j) is vaHd at £; < 40 TeV (with 7 = 1.62) and > 2 x 10^ TeV (with 7 = 2.02). The 
O (h/Xji) corrections are small at h < 500 g/cm^ and the derived solution will suffice for our purpose. 



E. Nuclear cascade at low and intermediate energies 



For the "low-energy part" of the nuclear cascade {Eq < 1 TeV/particle) , we adopt the relevant results of Refs. [t77|j78| 
obtained within a rather circumstantial nuclear-cascade model. The model includes the effects of strong scaling 
violation in hadron-nucleus and nucleus-nucleus collisions, ionization energy losses of charged particles, temperature 
gradient of the stratosphere, geomagnetic cutoffs and solar modulation of the primary spectrum. The computational 
results were verified considering a great body of data on the secondary nucleons, mesons, and muons in wide ranges of 
geographical latitudes and altitudes in the atmosphere. The model was also tested using low-energy data on contained 
events observed with several underground neutrino detectors. 

Since the key features of the model were discussed in several papers (see e.g. Ref. |6j and references therein), we 
shall not dwell upon the question here. Only one point needs to be made. The geomagnetic effects for the sea-level 
muon flux are sizable up to about 5 GeV. However, later on, we are going to deal with the muon data at high latitudes 
that are insensitive to the geomagnetic cutoff. The same all the more true of the solar modulation effects. 
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III. CONVENTIONAL MUON FLUX 



Our calculation of the muon production and propagation through the atmosphere is based on the standard continuos 
loss approximation. Our interest is in the muon flux at momenta p > 1 GeV/c. Thus the O (mf^/p-^) effects can be 
neglected. For simplicity, the nonisothermality of the atmosphere will be ignored in the formulas which follow (see 
Ref. for the corresponding corrections). 

Let 'D^{E, h, I?) be the differential energy spectrum of muons at depth h and zenith angle and Pfi{E) ~ —dE/dh ~ 
a^{E) + bfj,{E)E be the rate of the muon energy loss due to ionization (a^(_E)) and radiative and photonuclear 
interactions in the air [b^[E)E). The muon transport equation is 



_9_ 

dh 



Eh 



d 



V^iE, K^)^^ [P^iE)V^{E, h, d)] + Gl''^ [E, h, d) 



(3.1) 



with 




h.d\dx 



(3.2) 



Here E"{'d) — m^iJo seci^/r^ ~ 1.03 sec?? GeV is the muon critical energy, i?(M^2(3)) are the branching ratios for 
the 7r^2, K^2, and ii'^a decays, F^{x) is the muon spectral function for iiT^s decay, and 



2m: 



{m\ -ml+ ml) ± (mj^-ml + ml)' 



The explicit form of (x) is rather cumbersome but there is no need to write it out because the K^^ decay contribution 
to the muon flux does not exceed 2.5% | |79[ |. 
The solution to Eq. (3T) is given by 

V^{E,h,d) = 

Here 

W^,{E,h,h',d) = ^^-^ ^exp 



Wf,{E, h, h', ^)Gl'^ {£{E, h-h'), h', I?) dh'. 



/3^{£{E,h-h')) 



E'jm dh" 
y S{E,h-h") h" 



(3.3) 



and £{E, h) is the root of the integral equation 



dE 

ME) 



that is, the energy which a muon must have at the top of the atmosphere in order to reach depth h with energy E. 
As our analysis demonstrates, the weak (logarithmic) energy dependence of the functions and 6^ is only essential 
for near-horizontal directions and can be disregarded with an accuracy better than 3 % for the directions close to 
vertical. In this approximation 



£{E,h) = E 



exp(6^/i) 



and 



PA£iE,h)) 
ME) 



exp(&^/i). 



In numerical calculations we use — 2.0 MeV-cm^/g and 5^ — 3.5 x 10 ^ cm^/g. Eq. (3.3) can be simplified in the 
two particular cases. At i? » E^'^^{'d) the muon decay can be neglected, so 



W^{E,h,h',i))^exp[b^ {h~h')]. 
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At E <^ a^/h^ K, 0.57 TeV, the radiative and photonuclear energy loss are inessential and thus 

The combined results of our calculations for the vertical momentum spectrum of conventional muons at sea level 
can be summarized (with a 2 % accuracy) by the following fitting formula: 

(p, /i = 1030 g/cm^ I? = 0°) Cp-(''°+''i^°s''+^^i°s'f+'^^'°s'p\ (3.4) 

with parameters presented in Table || for a few momentum ranges [here p is the muon momentum in GeV/c and 
V^{p,h,^) = {j>lE)V^{E,h,d)]. 



TABLE II. Parameters of the fitting formula (p.4|) for the vertical energy spectrum of conventional muons at sea level. 



Momentum range (GeV/c) 


C (cm-^s-^sr-^GeV-^) 


70 


71 


72 


73 


1 ^ 9.2765 X 10^ 


2.950 X 10~^ 


0.3061 


1.2743 


-0.2630 


0.0252 


9.2765 X 10^ -^ 1.5878 x 10^ 


1.781 X 10"^ 


1.7910 


0.3040 








1.5878 X 10^ -f- 4.1625 x 10^ 


1.435 X 10^ 


3.6720 











> 4.1625 X 10^ 


10^ 


4 












Figure |^ compares our result for the vertical differential muon spectrum at sea level with the results of Volkova et 
al. |8^, Dar Butkevich et al. [|2j, Lipari [|3|, and Agrawal et a l. l84|. In this comparison, we used the fitting 
formulas from Rcfs. [^,0, and the corresponding tables from Refs. [^2|- ^4| 

In Table III, we show the ratio of each calculated spectrum from Refs.lSQ-M to ours for i? = 1, 10, . . . , 10^ GeV. 
The ratios are inside the wide range 0.75 1.48. In the momentum region from ~ 5 to 5 x lO'^ GeV/c, our result is in 
very good agreement with the recent Monte Carlo calculation by Agrawal et al. |8^: the discrepancy is less than 6 %. 
This is consistent with the uncertainties of both calculations caused by the uncertainties in the input parameters. At 
low energies (1 10 GeV), our calculation agrees closely with the fitting formulas by Volkova et al. pG] and Dar M|. 



TABLE III. The ratios of vertical differential spectra of conventional muons calculated by different workers to ours. 



Ref. 






E (GeV) 










1 


10 


10^ 


lo-" 


10* 


10" 


10" 


1 


3C 


1.010 


0.996 


1.135 


1.056 


1.189 


1.156 


1.483 


1 


~ 


1.001 


1.046 


0.958 


0.873 


1.023 


1.047 


1.405 




32 




1.015 


1.079 


0.909 


0.958 


0.902 


1.140 


I 


i3 


0.753 


0.820 


0.858 


0.823 


0.955 


0.923 


1.160 


1 


i4 


1.355 


0.992 


1.017 


0.938 









IV. CHARM PRODUCTION AND PROMPT MUONS 



The prompt muon and neutrino component of the cosmic ray flux originates from the decay of short-lived particles 
(mainly charmed hadrons , , , A+, . . . ) produced in interactions of cosmic rays with the atmosphere. For the 
last fifteen years, a lot of papers with calculations of prompt lepton production in the atmosphere have been published 
with very different outputs. Suffice it to say that the predicted energy at which the vertical sea-level PM flux becomes 
equal to that of muons from tt and K decays varies from ~ 20 TeV to ~ 10'^ TeV, depending on an adopted charm 
production model. Early works [^-p2[ were based on empirical ad hoc models for open-charm production and some 
extrapolations of the accessible (rather fragmentary) accelerator data to the orders-of-magnitude higher energies of 
the primary and secondary particles participant in cosmic-ray interactions. The successive works apply more advanced 
phenomenological approaches to the charm-production problem Pj5|,|93|-p6t or a set of parametrizations for the energy 
dependence of the inclusive cross sections those qualitatively describe the main features of some popular models for 
charm production 0j9^ . Let us glance off two recent approaches based on the perturbative QCD and the Dual Parton 
Model (DPM) |,|6[. 
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n 1 2 4 

10 10 10 10 10 10 

Muon Momentum ( GeV/c ) 



FIG. 1. Vertical differential momentum spectra of conventional muons at sea level calculated by Volkova et al. [gO|, Dar [ pl| , 
Butkevich et al. j82|, Lipari [Q, Agrawal et al. and in present work. 

Thunman et al. |^ apply a state-of-the-art model to simulate charm hadroproduction through pQCD processes. 
To leading order in the coupling constant, as, these are the gluon-gluon fusion [gg — s- cc) and the quark- antiquark 
annihilation {qq — > cc). The next-to-leading order, 0{as), contributions are taken into account by doubling the cross 
sections. To simulate the primary and cascade interactions, the authors use the well-accepted Monte Carlo code 
PYTHYA. Without going into details of their approach, we emphasize that the PM flux predicted by Thunman et al. 
is one of the lowest ones. It overcomes the vertical tt, i^-muon flux at energy of about 2 x 10'^ TeV and therefore is 
undistinguished in present-day ground-based and underground/ water muon experiments. 

In the paper by Battistoni et al. a new Monte Carlo calculation of the PM fraction in atmospheric showers 
was made using the DPMJET-II code based on the two-component DPM and interfaced to the shower code HEMAS. 
The calculation does not yield the absolute PM flux but, from the estimated prompt-to-conventional muon ratio, 
one can see a leastwise qualitative agreement with the result of Ref. In particular, according to the DPM, the 
prompt component overcomes the conventional one in the region of a thousand TeV (not reachable with the simulated 
statistics). 

In our previous works P j9^ , p^ , the two different phenomenological nonperturbative approaches to the charm- 
production problem have been applied, the Recombination Quark-Parton Model (RQPM) and the Quark-Gluon 
String Model (QGSM). In the present calculation, we use just these two models. For this reason, the most salient 
features of them will be outlined below in this Section. The RQPM will be discussed at greater length, considering 
that the QGSM is well accepted and covered adequately in the literature |Q (see also Ref. and Q for reviews). 
As an example of a calculation giving a particularly high PM flux, we will also sketch a semiempirical model put 
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forward by Volkova et al. ]92[ | . The comprehensive reviews of the current experimental status of the charm production 
problem can be found in Ref. |100|. 



A. Models for charm hadroproduction 

1. Recombination quark-parton model (RQPM) 

The RQPM is one of the models with "intrinsic charm" . The models of this class are based on the following key 
assumptions. 



(i) The projectile wave contains an intrinsic-charm Fock component (see Refs. |101,102]). As an example, Figure g 
shows the component \uudcc) generated by the virtual subprocess gg cc where the initial gluons couple to 
two (or more) valence quarks of the projectile. 



U 
U 



c 
c 
d 



FIG. 2. Intrinsic \uudcc) Fock component in the wave function of a projectile proton. 



(ii) The interaction of partons in the final state leads to a recombination (or coales cence) of the charmed quark with 
projectile fragments and to production of leading charmed hadrons [F03 Hl05| 



An indication in favor of these models was found in muon-nucleon scattering [ |106| . It was shown that there exists a 
visible excess of the charmed particle yield at xp ^0.15 and < 40 GeV^ over the model expectations based on the 
photon-gluon fusion and conventional QCD evolution. The upper bound for the probability to find an intrinsic-charm 
Fock component in the proton wave is abou t 0.6%. 

It has been shown by Brodsky et al. |102|] that the diagrams with intrinsic charm, in which a cc pair is coupled to 
more than one constituent of the projectile hadron, are suppressed by powers of M^-{1 — Xc) (here Mcc is the invariant 
mass of the pair and Xc is the fraction of hadron momentum carried by a parton), i.e. the relative contribution 
of the intrinsic-charm mechanism to the longitudinal momentum distribution of charmed hadrons is expected to be 
especially large in the fragmentation region of a projectile. In other words, intrinsic-charm models predict relatively 
hard inclusive spectra. At the same time, the total inclusive cross section can be rather large (it depends strongly on 
the assumptions about the charm structure function of the projectile hadron). These features cannot be obtained in 



perturbation theory (see e.g. Ref. |107] where a comparison of 600 GeV n emulsion data with the next-to- leading 
order pQCD predictions was made). 



In the RQPM, the process of hadronization occurs by means of recombination of quarks to hadrons [104|. It is 
assumed that only slow ( "wee" ) partons of colliding hadrons take part in the interaction and the distributions of fast 
partons do not change during the collision. Therefore the inclusive spectra of produced particles (those with small pT 
and with not too small xp) are entirely governed by quark distributions inside the projectile hadron. 

a. Charm production in hadron-nucleon collisions. Inclusive cross section for production of a meson M = qq in 
pp interaction is 

= I {^q^^^1j) ^pl^ {Xqi)Fpf {Xq-,Xq,Xg) Rm {Xq, Xq] X p) dXq^dx q.dXqdXq. (4.1) 
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Here qi and qj are the "wee" partons from protons pi and p2, respectively {p2 is the projectile), aij is the total cross 

section for qiqj interaction, -Fp™^ is the m-parton joint distribution ins ide the proton pk, and Rm is the function 
of recombination of the pair qq into meson M . The cross section ( [4.1[ ) is written for the fragmentation region of 
the projectile p2. Let us assume that the distribution of "wee" partons is universal and does not correlate with the 
distribution of fast partons. Then 



P2 



^li ) ^P2 (^9' ^9) 



Considering that 



yields 



E 



Qj ' 



pp^MX 



F'^^ (Xg, Xq) i?M {Xq, Xq-, X p) dXqdXq. 



dXF 

In a similar spirit one can derive the inclusive cross section for the generic reaction iN — )■ fX: 

da,N-.fx ^^tot(^) J F,{{xk})Rf{{xk};xF)l[dxk. 

k 



dxF 



(4.2) 



Here x^ is the fraction of the projectile momentum which belongs to the parton g^, ({xfc}) is the two- or three-quark 
distribution in the projectile hadron i and Rf {{xk\] xf) is the function of recombination of two or three quarks into 
hadrons. 

It would appear reasonable that far away from the threshold of open-charm production, the parton distributions and 
recombination functions do not depend on the projectile particle energy. Then the s-dependence of the daiN-,fx/dxF 
is determined by the energy dependence of the total cross section for the iN interaction, tT'^(s), and therefore the 
scaling violation is fairly small. As in the case of light particle production, we use for the cr*]^(s) the model of elastic 
amplitude from Ref. ||7J which predicts that the total cross section grows as In s at the asymptotic energies. 

We assume that the c- quar k sea in a hadron is essentially nonperturbative and it is characterized by a flat momentum 
spectrum (see e.g. Ref. [101|). According to the parton conception, in the infinite momentum frame, the lifetime of 



fluctuations containing heavy quarks is very large; the flatness of heavy-quark spectra follows from a simple picture of 
a hadron as an aggregate of partons with approximately equal velocities and from calculations of structure functions 
for strongly coupled states. 

To calcu late the two- and three-quark distributions, Fi {{xk}), we use the statistical approach by Kuti and Weis- 



skopf |l08[ . The functions Fi {{xk}) are constructed through "uncorrelated" parton distributions fk'^^{x) and fl°'^{x) 
for valence and sea quarks (fc = u, d, s, c) and through the correlation functions G^{l~x). For example, the two-particle 
distribution of u and c quarks in a proton is of the form 

(^^,2;^) = [2GP(1 -x^- x^)f:^\x^) + Glil - x„ - x^)fnx^)] frix^)- 

The three-particle distribution of u, d, and c in a proton is 



(x„, X,, X,) = [2GI, (1 - 5] Xq) /r'(x„) + (1 - 5] Xq) frix^)] fT\x,)fr{xc) 
2Gi (1 - ^ Xq) f:^Hxu) + (1 - ^ Xq) frix^)] frixd)fri^c), 



Xu ^" Xfl Xq. 



Both the uncorrelated distributions and correlation functions for light quarks and gluons in a proton and pion were 
calculated by Takasugi [109| in the framework of the statistical model using all appropriate accelerator data. It can be 
shown that the correlation functions are little affected by introducing the sea of charmed quarks and hence we will use 
the results of Ref. |109{ without any modifications. In so doing and using Eq. (4.2), the uncorrelated c-distributions, 
fc°^{x), could be basically extracted from the data on charm production. In fact the realization of this program is 
somewhat limited because Eq. ( |4.2| ) only holds at asymptotic energies (far away from the charm-production threshold) 
and besides, the available accelerator data at high energies cover a narrow range, 0.1 < a;_F < 0.9. Within this range, 
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the best fit of the ISR data o n production in pp interactions |110| and the EMC data on charm production in 
deep- inelastic muon scattering |l06t| is achieved with the following simple parametrizations [ 104 1 : 



5.5 X 10"^a;""-^(l - x)-^-^^ for proton, 

\-0.85 



7.7 X 10 



-3^-1 



(1 



for pion. 



In our calculations, we do not make distinctions between pseudoscalar and vector charmed mesons of an identical 
quark composition at production. So, by a Z? meson production cross section is meant the overall cross section for 
production of D and D* mesons. 



For the recombination functions of quarks into D and Ac we use the formula derived by Hwa in his valon model [111 



Rd{xi,X2]x) 
RkSxi,X2,xz\x) 



B{a,b) 



^(-^.^ - a;), 



B{a,b)B{a,a + b) 



mi 



I 6{xi 
x / 



+ X2+XZ- x). 



Here B{a, b) is the beta-function, a and b are the constants defined by the form of t he v alon distributions. Regarding 
the valons as constituent quarks bound nonrelativistically in a bag, it can be shown |111| that their average momenta, 
(xi), are proportional to their masses, rhi. Then, considering the two- valon distribution in a D'^-meson, we have 



a 
b 



rrh. 



Below, we adopt a = 1 in all numerical calculations. 

h. Nuclear effects. In order to take the nuclear effects into account, we use the additive quark model 111 



Let us 

assume that passing over the target nucleus, a valence quark of the projectile behaves as a free particle between 
its collisions with nuclei. If at a collision with a nucleus the quark loses the bulk of its momentum, that quark may be 
thought of as captured by the target and its contribution to the production (through the recombination) of hadrons 
with large xp can be neglected. On the contrary, the quark which escape collisions can hadronize by recombining with 
slow quark(s) as described above. Because our prime interest is in the high-energy range and in the fragmentation 
region of projectile particles, one can neglect the interaction of secondary hadrons with the target nucleus. Indeed, the 
time of generation of hadrons is proportional to their momenta and fast hadrons are produced outside the nucleus. In 
line with these assumptions, the invariant cross section for inclusive production of hadrons in hadron-nucleus collisions 
is expressed in terms of the "recombination" hadron-nucleus cross sections and the probabilities for capturing valence 
quarks by the target nucleus . Us ing standard "nuclear optics" techniques |113] and the additive-quark- model relations 
for the total cross sections [112], 2app ~ Sav^ — 2a qp [q — u,u,d,d), one can derive the following formulas for the 



inclusive charm-production cross sections [104[ 



dxF 



dxF 

d<JpA^DOX 



dxi 



do-pA^D-X _ Q / '^PA ~ '^1^ \ pp^D-X 



dxi 



2 l' UtvA - CFqA \ "''^pp^D-X 



Jpp 



dxj 



dxp 



~pp-*D°X 

dxi 



da 



dxp 



fpA + (^ttA — "^f^qA \ ' ^pp^WX _|_ 3 /' ^ttA — <^qA \ pp^T5° X 



da 



pA^htx 

dxp 



^ pp 

2 / <7pA ~ CTttA 

3 / apA - aqA 
2 I 



I ff 
J dxp 



Jpp 



da^ 



■ 



a 



pp 



dxp 

1 *- 

J dxi 



Jpp 



da 



dxi 



da^ 



da^+A-^DX _ c) ( ^tA — <^qA\ da^yp^DX 



dxi 



dxp 

pp^A+X , 2 / <^ttA - <^qA \ "'^pp^KtX 

app J dxp 
- n± nO 7T0^ 



dxi 



[D = D^,D°,D°). 



Here c?ct|p_^ /dxp is the contribution to the iN cross section from a quark diagram with a final hadron / that contains 
the leading valence (v) or sea (s) quarks indexed in the brackets. To sufficient accuracy, the total cross sections in 
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the foregoing equations are assumed to be energy- independent. In our numerical evaluations, we set aiA = ct ^a for 
the hadron-nucleus cross sections (see Section ||) and cr^p = 13.0 mb for the quark-proton cross section [113|. The 
numerical results are represented in the traditional form, 

For the reactions pA -> D+X, pA D°X, and nA DX {D ^ D'^,D",D^), a = 0.765, independently of xp- It 
should be pointed out that accelerator data at low energies show a higher value of a. For example, in the WA82 



experiment 114 (a 340 GeV tt" beam) the value a = 0.92 ± 0.06 was obtained for D mesons with (xp) — 0.24. 
However, it seems plausible that this is a reflection of the "near-threshold effect" and the a will decrease with a rise 
of the projectile energy. In any event, the non-perturbative effects should become more important as ^/s and xp 
increase and t herefore the shadowing is expected to become more essential at higher center-of-mass energies and at 



large xp |]115| 



For the rest reactions and within the range 0.10 < xp < 0.95, the functions a{xp) may be parametrized as follows: 

%A~.D«xi^) = 0-754 - 0.034a; - O.OOSa;^ -I- 0.0202:^ 
apA-^D-x{x) = 0.769 - 0.158a; + {).212x^ - 0.174x^ 
%A~,Atx^^) = 0-''80 - 0.367a; + 0.6722;^ - 0.456a;^ 
These results do not contradict the acce lerat or data even at very low energies, although the data are rather uncertain 



yet. For example, the BIS-2 experiment |11£] (a 37.5-70 GeV neutron beam) gives (a) = 0.73±0.23for D° production. 

As discussed above, we assume that the captured quarks take no part in the recombination. This leads to a 
small underestimation of the cross sections, because some portion of wounded quarks actually will recombine. Let 
us estimate the upper limit of the a assuming that all the valence quarks of the projectile can recombine. This 
assumption yields 

dfTiA^fX _ ( OiA \ d(7iN~*fX 



dxp ycTiN J dxp 

and thus a < 0.85 for ttA DX and a < 0.79 for pA D {Ac)X. This estimate demonstrates that the uncertainty 
in the A-dependence within our simplified approach does not exceed '--^ 15% for the air nuclei. 

c. Z-factors. Ow ing to the mentioned small scaling violation, the fractional moments Zfi calculated with the 
RQPM from Eq. ( |2.5| ) are energy dependent. They can be approximated with an accuracy of (2-3) % by the following 
expression: 

Zf,{^,E) = Zf.,{j,E,)(^-^^^\ (4.3) 
where E is the energy of secondary particle f {f — D^, D", D*^, A+), E^ and arc the constants dependent on the 



TABLE IV. Parameters Zfi['y,E-y) of fitting formula (4.5) for the fractional moments Zfi{-f,E) calculated with the RQPM 
for the two values of 7. 









/ 






i 


D+ 


D- 




W 


A+ 








7 = 1.62 






P 


4.6 X 10"* 


6.5 X 10"* 


3.8 X 10"* 


6.9 X 10"* 


4.9 X 10"* 




1.3 X 10"^ 


9.0 X 10"* 


9.0 X 10"* 


1.3 X 10"^ 


6.0 X 10"* 








7 = 2.02 






P 


5.4 X 10~* 


7.9 X 10"* 


4.5 X 10"* 


8.6 X 10"* 


6.2 X 10"* 




1.8 X 10~^ 


1.2 X 10"^ 


1.2 X 10"^ 


1.8 X 10"^ 


7.9 X 10"* 



primary cosmic-ray spectrum. In particular, 

E-y = 1( 

Ey = W GeV, ^-y = 0.076 for 7 = 2.02 



Ey = 10^ GeV, ^-y = 0.096 for 7 = 1.62, 
-16 
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Parameters Zfi^j, E^) for i = p, tt"*" are presented in Table [V. For i = n,?: one can use the relations 



— ^7=70 , 
D p' 



- Zo-p, ^A+n - -^A+p' 



which follow from considerations of the isotopic symmetry. 



A+7r+' 



2. Quark-gluon string model ( QGSM) 



The QGSM |98| is a non-perturbative approach to the description of hadron collisions. It is based on the topo- 
logical 1/Nf expansion of QCD diagrams for elastic scattering |117| (associated with the multiple pomeron exchange 
expansion) and the string model of hadrons and hadronic interactions. The particles are produced in this model by 
breaking the strings connecting the incident hadron's constituents (quarks and diquarks). 

The QGSM is considered to be one of the most satisfactory of the tools available to represent open-charm production. 
It describes a great body of data on hadronic interactions at all available energies. However, the model is not free from 
difficulties. For instance, the QGSM predicts clear-cut flavor correlations. In particular, there must be preferential 
production of D'^ mesons in pp collisions ("favored fragmentation") owing to (u — u d) co mposition of the proton and 
(cm) composition of D*^ (Figure^). This prediction is not supported by experiment |118|, although this disagreement 
can be caused in part by bad flavor identification in the experiment (see Ref. M for a discussion). 



U 



u 



ud 



cc 




ud 



cc 




u 



uu 



cc 



dd 



u 




(a) (b) (c) 

FIG. 3. Fragmentation of quark chains into D mesons in the QGSM: (a,b) favored fragmentation into D''; (c) unfavoured 
fragmentation into D~ and D". 

To calculate the inclusive cross sections one must know the distribution functions of the dressed quarks (constituents) 
of the colliding hadrons and the fragmentation functions of these constituents into charmed particles. These functions 



can be approximately determined by the use of Regge model arguments [119|, in terms of intercepts an —oim ~ 0.5, 
of known Regge poles and the intercept of the cc Regge trajectory, a^, on which there is no direct experimental 
information. Hence is a free parameter of the model. It governs, in particular, the steepness of the inclusive 
spectra of charmed particles. If the cc trajectories are linear (as it is in the case of light quarks and generally in the 
string models of hadrons), the intercept of the 7p trajectory is fairly large (~ —2.2) and the longitudinal momentum 
distributions of charmed hadrons are rather steep. A complete list of the distribution and fragmentation functions as 
well as the values of their various parameters are given in Ref. ||98|| . 

Our calculations of the inclusive cross sections within the framework of the QGSM have been done without attempts 
to optimize the set of parameters of the model. In particular, we do not include the intrinsic charm component as 
it was suggested recently . Below, we are dealing with a qualitative analysis of the QGSM prediction for charm 
production at cosmic-ray energies rather than with a close examination of the model. For this reason, in evaluating 
the nuclear effect within the QGSM, we adopt a = 0.72 for all processes under consideration. This simplification can 
lead to a small (< 15%) error in the Z-factors, compared to the exact calculation within the additive quark model. 
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The energy dependence of the factors Zfjj^, E) calculated with the QGSM is somewhat different as compared with 
the RQPM prediction. The parametrization ( [4.3|) is vahd for the QGSM only at very high energies (> 10"^ TeV) and 
the parameters £,j are in general different for different reactions iA fX. The parameters Zfi{'-f,E~^) and for 
i = p,n, 7r+ and 7r~ are presented in Table ^ at 7 = 2.02 (above the knee energy region). The energy dependence of 
the Z-factors at £; < 10^ TeV can be found in Rcf. 101. 



TABLE V. Parameters Zfi{^,E^) and ^. 
calculated with the QGSM for 7 = 2.02 at E > 10^ TeV 



(in parentheses) of fitting formula (4.3) for the fractional moments Zfi{'y,E) 

,3 









f 






i 


D+ 


D- 




D° 




P 
n 

7r+ 
7r~ 


6.5 X 10-" 

(0.050) 
7.1 X 10"^ 

(0.050) 
5.5 X 10"'' 

(0.041) 
1.4 X 10"* 

(0.041) 


9.9 X 10"'' 
(0.046) 

1.9 X 10"" 
(0.045) 

1.4 X 10"'' 
(0.048) 

5.5 X 10"* 
(0.048) 


7.1 X 10"^ 
(0.050) 

6.5 X lO"'^ 
(0.050) 

1.4 X 10"" 
(0.048) 

5.5 X 10"* 
(0.048) 


2.1 X 10"" 
(0.044) 

1.2 X 10"'' 
(0.045) 

5.5 X 10"* 
(0.041) 

1.4 X 10"* 
(0.041) 


9.5 X 10"" 

(0.041) 
9.5 X 10"* 

(0.041) 
1.5 X 10"^ 

(0.035) 
1.5 X lO"'^ 

(0.035) 



3. Semiempirical model (VFGS) 

The model of Volkova et al. (let us call it the VFGS model) is a typical example of an approach which proceeds 
from a parametrization of available accelerator data for inclusive spectra of charmed particles together with some 
additional assumptions to extrapolate the parametrization to the kinematic regions, where the data on the inclusive 
charm production cross sections are absent. 

Volkova et al. make use a very steep inclusive spectrum of produced D-mesons (oc (1 — xd)'^/xd, where xo is 
the ratio of the D-meson energy to the nucleon energy in the lab. frame) with a sharp cut-off in the central region 
{da/dxD = at xd < 0.05). In spite of such cut-off the integral j [da / dxD)dxD was normalized to the total DD cross 
section, f7^^(i?jv). Considering the accelerator data at > 1 TeV together with some implications of the QGSM, 
it has been adopted that 

DD,j^ , r 0.48(log£:A, - 3.075) mb for 1 TeV < En < 500 TeV, 
CTpp {En) - I 26 mb for En > 500 TeV. 

A consequence of this assumption is a relatively strong scaling violation in the fragmentation region. 

The VFGS model predicts comparatively large PM flux (see below) since, owing to the cut-off, all produced particles 
are in the fragmentation region of a projectile (i.e. there is no the central part of the inclusive spectrum). It was also 
assumed that (independently of xp) ct = 1 and 2/3 for reactions with D mesons and hyperons in the final state, 
respectively. 

The approach of Ref. [ p2[ includes some other assumptions which also tend to increase the PM fraction in comparison 
with our result. The most important ones are concerned with the primary spectrum, semileptonic decays of charmed 
particles and with certain elements of the nuclear-cascade model. A more detailed comparison of the approach under 
consideration against the RQPM and the QGSM, in connection with the PM problem, has been done in Ref. |Q. 



B. Prompt muon flux at sea level 

1. Interactions and decay of charmed particles 
As we neglect the production of nucleons, pions, and kaons by c harm ed particles and charm regeneration, the 



transport equations for D and Aj. spectra are identical in form to Eq. (2.7) for kaons. Notice that the PM flux weakly 
depends on the specific values of the inelastic cross sections for D and Ac up to about 10* TeV of muon energy, due 
to very short lifetimes of these particles. Thus a rough estimation of ctJ^^' and tT™±'^ will suffice for our purposes. We 



use the same formula ( ^.4[ ) as for the light hadrons with ct^,^ = 100 mb {D = , D^, D^) and <^^±^ = 200 mb. 
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Calculation of the PM flux can be performed in almost perfect analogy to the conventional muon fluxe with the 
only one essential difference: the PM generation function includes a rich variety of multiparticle semileptonic decay 
modes. Thus the inclusive approach is best suited to the problem. The corresponding muon generation function may 
be written as 



J2 Bit -> ^ii^X) 

i=D± .D° .D« ,Aa 



hE 



(4.4) 



Here F^(x) is the normalized spectrum of muons in the inclusive decay i —t ^v^X [x = E/Ei) and 

= 2m^ (mf + ml - sx) ± \J (mf + m2 - sxf - 4TO2y 

with Sx the minimal invariant mass square for the hadron system X. The other designations are completely similar 
to the ones previously used. 

To simplify matters we consider the inclusive decay i fivX as a 3-particle one. We assume the simplest form of 
matrix elements according to Ref. [120|. The form factors involved (one for D — > ijlv^X and three for Ac fjiv^X) 
are replaced with their averaged values. In so doing the mass square of the "X-particle" , s^f i may be fitted in such a 
way as to correlate the calculated and experimental values for the differential and total decay rates. Omitting rather 
tedious details of the calculation, we present the final formulas for the muon spectral functions F^(x) and Fj^ (x) in 
Appendix ^. 



2. Parametrization of the calculated PM flux 

In the energy region 5 TeV ^ E < 5 x lO'^ TeV the differential spectra of PM in the vertical direction at sea level, 
Vf^'iE), calculated in Ref. ^ with the RQPM and the QGSM, can be approximated by 



(E,h = 10S0g/cm\^ = 0°) =C' ( ^ 



7-1' 



(4.5) 



Here 



C" = 4.53 X 10"^* cm"2s"^sr"^GcV"\ 7' = 2.96, 
C = 1.09 X 10"^** cm"2s"isr"iGeV"\ 7' = 3.02, 



a 
a 



0.152 (RQPM), 
0.165 (QGSM), 



and Eb = 10^ GeV in both cases. Eq. ( fl.Sj ) fits the numerical results with accuracy better than 4 %. With the same 
accuracy it is also valid for zenith angles -d < 80° in the energy interval (10 10^) TeV, i.e. with in t he "region of 
isotropy" of the PM flux (see Ref. H for more details). Beyond the interval (5-f 5 x 10'^) TeV, Eq. ( |4.5| ) can be used 
as an extrapolation of our result which would suffice for calculating the muon DIR. 



It is interesting to note that the RQPM and QGSM predict very different values for the muon charge ratio [121 
The energy dependencies of the charge ratios may be approximated by 



V 



pr 



fJ- 



0.864 - 0.006 log^ {E/En) for RQPM, 



1.250 + 0.008 {E/Er 



,0.73 



for QGSM, 



with Efj = 10 TeV. These approximations are valid in the energy range 3 10'^ TeV at all zenith angles with an 
accuracy bet ter t han 2%. 

From Eq. (4.5) we find the following expression for the integral PM spectrum: 



JP' [E,h^ 1030 g/cui , i9 = 0' 



C'Eb 



(7'-l)(l-a) 



El 
E 



A comparison of our calculation of the PM flux with the results of other authors can be found in Refs. (see 
also M). 
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According to Ref . [g^ , the differential and integral PM spectra calculated in the VFGS model can be approximated 
(at all zenith angles) by 

{E, h = 1030 g/cm^ I?) = 2.92 x IQ-'^E-^ *^ cm"2s"^sr"iGeV~\ 
IP'- [E, h = 1030 g/cm^ z9) = 1.97 x IQ-'^E^^ *^ cm^^s^^sr"^ 

{E in GeV.) This approximation holds true to about 10"^ TeV. 

V. CALCULATED SEA-LEVEL MUON SPECTRA VS EXPERIMENT 

Comparison of the calculated differential and integral muon spectra with direct data from spectrometers and 
indirect data extracted from underground measurements is shown in Figures^ (a,b) and ^ (a,b). The ground-based 
measurements can be classified as absolute and non-absolute (normalized). In line with this arrangement we present 
here the following three groups of experiments. 

• Absolute ground-based measurements 

with MARS apparatus in Durham (Aurela et al. 0, Ayre et al. 12 1); Nottingham spectrograph 
(Baber et al. ||], Rastin jlj]); spectrometer near College Station, Texas (Bateman et al. [^); Kiel 
spectrographs (AUkofcr et al. ||l^), MASS apparatus at Prince Albert, Saskatchewan (De Pascale et 
al. [^); EAS-TOP array at Campo Imperatore, Gran Sasso (Aglietta et al. jT7|]). 

• Non-absolute ground-based measurements 

with Durgapur spectrograph (Nandi and Sinha ]ll[ , the data were normalized to the Nottingham 
spectrum |lj] at p = 20 GeV/c); Durham spectrograph MARS (Thompson et al. ||l^, the data were 
normalized to the previous MARS results ||l2| at 261 GeV/c); L3 detector at CERN, (Bruscoli and 
Pieri |l^, the absolute intensity in the momentum range 40-70 GeV/c and its error were taken from 
the Kiel result 0). 

• Indirect data 

from several detectors in the Kolar Gold Fields (Ito Miyake et al. [Q, Adarkar et al. psf); 
unimodular scintillation detector "Collapse" of the Institute for Nuclear Research (INR) at the Arty- 
omovsk Scientific Station (Khalchukov et al. ) ; Baksan underground scintillation telescope of INR 
situated in North Caucasus (Andreyev et al. , Bakatanov et al. |64j ) ; X-ray emulsion chambers 

of Moscow State University situated in the Moscow metro (Zatsepin et al. |^^); proton decay detec- 
tor Frejus under the Alps (Rhode ||6^), detector MACRO at the Gran Sasso National Laboratory 
(Ambrosio et al. Q). 

The marked curves in Figures ^ and ^ refer to the differential and integral muon spectra, respectively, calculated 
without the PM contribution ("tt, _ftr"-muons) and with the PM contribution according to the three charm production 
models (QGSM, RQPM, and VFGS) under consideration. As seen from the Figures, the PM contribution to the 
sea-level muon flux calculated with the QGSM is very small: up to p = 100 TeV/c it does not exceed 16% for the 
differential spectrum and 22 % for the integral spectrum. 

Unfortunately, it is difficult to extract some quantitative assessment for the validity of our nuclear cascade model 
from the presented set of data even at p < 1 TeV/c. As is seen from Figures ^(a) and H(a), a wide disagreement 
between the results of different experiments takes place despite the fact that the quoted errors are relatively small in 
the majority of the experiments. It indicates the existence of significant systematic errors in some experiments which 
may be as much as (30-35) % at momenta 10 to 1000 GeV/c. 

It should be noted in this connection that only statistical errors are indicated in the data points of the MASS 
experiment. According to Ref. [^, the systematic errors in the MASS experiment may be as much as 15% at 
p > 40 GeV/c. The systematics in the non-absolute measurements is, as a general rule, unknown. For example, 
no attempt was made to estimate the systematic errors in the CERN L3 experiment p6| . In our opinion, the L3 
spectrum was underestimated owing to incompletely correct normalization. 

At p < 2TeV/c our prediction, regardless of the charm production model, is in very good agreement with the 
Nottingham direct and absolute measurements |p^. 
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FIG. 4. Vertical difFerential momentum s pec trum of muons at sea level. The direct data are taken from Refs. |8|-|l6|| and 
indirect (underground) data are from Refs. ]44[0,||j64|-|6) . The shaded areas are for the MACRO fit @. The solid curves 
represent the results of this work for the conventional (tt, K) differential muon spectrum and for the n, K muon spectrum plus 
the PM contribution calculated according to QGSM, RQPM, and VFGS. 
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FIG. 5. Vertical integral momentum spectrum of muons at sea level. The direct data are taken from Refs. [[^,p| |ll| , p^Jl4| , |l5| , ^ 
and indirect (underground) data are from Refs. p7| , ^8|j6o| , |6^ . The solid curves represent the results of this work for the 
conventional (vr, K) integral muon spectrum and for the tt, K muon spectrum plus the PM contribution calculated according 
to QGSM, RQPM, and VFGS. 
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At energies above a few TeV we only have indirect data at our dispos al [ 122| a nd the uncertainties (both statistical 
and systematic) are vastly greater here. The data of Refs. ||37| , t38| , [44|60| 6»j|j65|] have been deduced from the muon 
DIR measured i n di fferent rocks (Baksan, Kolar, Alpine, Gran Sasso). We will dwell on the initial underground 



data in Section VIl . Here, it should be pointed out that in all underground experiments, among the systematic 
uncertainties related to inhomogeneities in density and chemical composition of the matter overburden, topographical 
map resolution, muon range-energy relation, muon range fluctuations, effective differential aperture of the array, etc., 
another uncertainty is essential. It results from the necessity to assign some model for the energy spectrum and 
zenith-angle distribution of muons at sea level which are functions of the PM fraction in the muon flux or, to be more 
specific, the ratio X of prompt muon spectrum to the tt + K production one. Hence one is forced to assume some 
value of the ratio X (as a function of energy) when reconstructing the vertical muon spectrum on surface. But the 
greater the adopted value of X , the harder the resultant spectrum. For this reason alone the conversion procedure is 
fairly ambiguous. 

As an illustration we consider the KGF results. The KGF muon spectrum in the energy range (200 7500) GeV 
was deduced Q using the underground data from Rcf. a nd assuming X = 0, what is quite reasonable for this 
range. But the data at higher energies (see also Rcf. ||63|l) demand a nonzero X. To estimate the ratio X, the 
authors have assumed a pion production spectrum of the form F{Et^) oc E^^ and a K/tt ratio of 0.15. The X ratio 
was assumed to be a constant. Then a analysis indicated that with 7 = 2.7 for muon energy of 8 to 250 TeV, there 
is PM production at the level of X = (9 ± 2) x 10~*. In Figure ||(b), we show this result (the corresponding data 
points are represented by diamonds) together with the spectrum deduced on the assumption that X — Q (the data 
points are represented by symbols x). As would be expected, the spectrum reconstructed with X = is softer. It 
is not difficult to understand that the final result is subject to variation also in response to variation of the adopted 



Kj-K ratio and 7 [123|. It should also be recognized that the real spectra of muons and mesons are far short of being 
power-law ones. 

Let us touch briefly on some essential points of the rest of the underground data presented in Figures (b) and 

1(b). _ 

In the Baksan experiment X = (1.5 ± 0.5) x 10 was found as the best fit of the calculated total intensity 
of conventional and prompt muons to the experimental data, assuming a power-law primary nucleon flux with the 
spectral index 7Ar = 1.65. 

In Ref. the complete data set of downgoing muons recorded with the Frejus detector has been reanalyzed. 
However in this analysis, the sea-level spectrum was derived using in essence the continuous loss approximation with 
some effective and energy-independent energy loss coefficients. The muon range ffuctuations are discussed in Ref. |6|] 
exclusively to estimate the uncertainty of the analysis. But it is a matter of common knowledge that, on calculation of 
the muon DIR, the continuous loss approximation results in downward bias and the corresponding error increases fast 



with depth [124-127|. It is our opinion that the muon spectrum obtained in Ref. [|65[ was significantly overestimated 
while the systematic errors were underestimated for > 10 TeV in consequence of the oversimplified analysis. 
The MACRO fit [H presented in Figs. 0(a,b) by shaded areas has the following form: 



Vr-° (e.H^ 10=g/c■n^.) . C„ ' + , (5.) 

^ ^ \ 115 GeV 850 GeV / 

with Co = (0.26 ± 0.01) cm^^s^^sr^^GeV^^ and 7^ — 2.78 ± 0.01. The quoted errors are due to statistics and 
the topographical map resolution. According to Ref. Q, the overall systematic error resulting from rock density 
uncertainties and hard energy loss cross sections is about 5% in Co and, what is much more important, 3% in 7^. 
But a 3 % variation in 7^ corresponds to uncertainties of 47 %, 78 % and more than 100 % in the surface muon flux at 
energies of 10^, 10'^ and 10* GeV, respectively. Therefore, the result of MACRO is greatly uncertain and pro forma 
it is not in contradiction with all the rest of indirect data shown in Figure ^. 

The results of the rest of the underground experiments, were obtained with quite different methods. The experiment 
with the Artyomovsk 100-ton installation "Collapse" |g2[ (situated in a salt mine at the depth of 570 m w.e.) detects 
the energy release of the showers produced by cosmic-ray muons in the salt and scintillator (Cioi?22)- In the Baksan 
"calorimetric" experiment ||6^ , the integral muon intensity at the position of the scintillation telescope (8.5 km w.e.) 
was evaluated from the spectrum of electromagnetic cascades generated by muons in the telescope. To find the muon 
intensity at the surface, the authors used a conversion procedure similar to that which was used in Refs. Due 
to a 10 % error in the calibration of the energy evolution in the detector, the systematic error in the determination of 
the absolute muon intensity can reach 25 % in this experiment. One might expect a supplement systematic uncertainty 
due to the conversion procedure. Comparing with the results of other experiments, the authors moved up their data 
by 12%. We use the same normafization in Figs. |(b) and |(b). The data of Moscow State University (MSU) |6^] 
were extracted from a multidimensional analysis of the measured energy and angular distributions of electron-photon 
cascades generated by muons in X-ray emulsion chambers. However, the output of this method is also very sensitive 
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to the adopted models for the primary spectrum and charm production. According to Refs. ^^l; the estimated 
primary spectrum index is = 1.64 ± 0.03 at nucleon energies (20 400) TeV and the best-fit X ratio changes from 
(2.6 ± 0.8) X 10-3 at E^5 TeV to (3.3 ± 1.0) x 10"^ at £; = 40 TeV. 

Up to about 5 TeV/c, our prediction for the differential and integral spectra (irrespective of a charm production 
model) does not contradict to the results of the Artyomovsk detector and the X-ray emulsion chambers of MSU. 
Below a few TeV the predicted spectrum agrees well with the data from Baksan, Frejus and MACRO extracted from 
the muon DIR, as well as with the Baksan data obtained from the spectrum of electromagnetic cascades JM. 

The region from 5 to 15-20 TeV/c is rather oracular: the data of KGF , Artyomovsk Baksan [|64[ and one 
data point of MSU |Q show a broad dip in the differential and/or integral spectra, whereas the rest of the data 
indicates some flattening or even a bulge 

Above ~ 20 TeV/c, the data of Baksan |37[, MSU and Frejus clearly indicate a significant flattening of 
the muon spectrum. Neither the QGSM nor the RQPM can explain this effect; even the maximum VFGS flux is not 
sufficiently larg e to this end, although the VFGS flux is not in contradiction with these data. It will be demonstrated 
in Section VII that this flattening is not confirmed by the body of direct underground data, while the late result of 
KGF seems to be (somewhat) more credible. It is also of interest that, irrespective of a charm production model, 
our prediction for the horizontal muon spectrum is in agreement with the corresponding MSU data | |6^ up to about 
40 TeV/c. The KGF spectrum (H] obtained at X = (9 ± 2) x 10""^ is in qualitative agreement with the RQPM 
prediction. Apparently the inconsistency of the data from different experiments gives no way of deducing a definite 
conclusion about the PM fraction in the sea-level muon flux. 



VI. MUON PROPAGATION THROUGH MATTER 



To calculate the muon depth-intensity relation we apply the semianalytical method proposed in Ref. | 127 |. The 
method allows us to avoid any simplifying assumptions about the scale invariance of the cross sections for radiative 
(direct e^e~ pair production, bremsstrahlung) and photonuclear interactions of muons with matter, and to take into 
account the real non-power-law behavior of the muon boundary spectrum. The solution to the transport equation for 
the differential muon intensity, 'D^{E, h), is constructed by iterations, starting from an initial approximation with the 
correct high-energy asymptotic behavior. Let us sketch the basic ideas and formulas. 

The equation describing the high-energy muons propagation through a homogeneous medium may be written 

o O /> 1 

—V^{E,h)^ — [B{E)V^{E,h)]^ J [il~v)-^^kiv,E,)Vf,iE,,h)-<i>k{v,E)V^{E,h)]dv, (6.1) 

with the boundary condition T>^[E, 0) = 'Do{E). Here 'Dq{E) is the ground-level muon spectrum, B is the rate of the 
muon energy loss those are treated as continuous. In the present calculation, B includes the ionization energy loss 
and the part of the loss due to e+e" pair production with v < vq — 2 x 10^**, where v is the fraction of the energy lost 



by muon (see Appendi x M and Ref. [128|). However the method is independent of the specific choice of B{E). The 
right-hand side of Eq. (|6.l[ ) describes the "discrete" muon energy loss resulting from direct e^e~ pair production with 
V > vo {k — p), bremsstrahlung (fc = b), and inelastic nuclear scattering (fc — n). The corresponding macroscopic 
cross sections, ^k{v,E), are defined by 



dv " dE' 



where A^'o is the number of atoms per 1 g of the matter and E {E') is the initial (final) muon energy. It is implied that 
the differential cross sections are averaged over the atomic number and weight of the target nuclei and dak {v, E) j dv = 
outside the ranges < v^™[E) < v'^'^^{E) < 1 allowed by kinematics. Lastly, Ey = E/(l — v). The summary of the 
explicit formulas for the cross sections used in our ca lcula tion is presented in Appendix 
Let us seek the solution of the transport equation ( |6.lD in the form 

V^{E, h) = VoiSgiE, h)) exp [~JC{E, h)] [1 + 5(S, h)]. (6.2) 

The functions involved are defined by the following chain of equations: 

^.j, f-^^"^ aE')~C{E')giE')~B'{E') 
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I. Jo I, Jo 



Evm 



VojEy) 
{l-v)Vo{E)- 



[As is easy to see, ( = "/+! and r/(u, E) — {1 — v)^ in the special case of a power-law boundary spectrum, VoiE) cx 
^-(7+1) ] xhe function £g{E, h) is the only root of the equation 



dE' 



E B{E') + q{E')E' 



(6.3) 



it can be treated as the effective energy, which a muon must have at the boundary of the medium in order to reach 
the depth h having energy E with a nonzero probability. Lastly, the function S{E, h) satisfies the equation 



with 



n{E,E,,h) 



d{E,h)^y] f ^k{v,E,){n{E,E,,h)[l + S{E,,h)] 
- [1 + uj{E, h)v][l + 5{E, h)]} ri{v, E)dv, 

Vo{E) Voi£g{E,,h)) 



Vq{E,) Voi£s{E,h)) 



exp [IC{E,h) ~ IC{Ey,h)] 



(6.4) 



uj{E,h) 



Q{E)~Q{8,{E, h)) 
B{E) + g{E)E ' 



Q{E) = \i{E) - B\E)\E + aE)B{E). 



Clearly S{E, 0) ~ 0. We shall seek the solution to Eq. (6.4) using an iteration procedure. It is based on the following 
consideration. 

Let us suppose that the functions $fc(w,i?) and CiE) become energy-independent as E 00. If so, it is a matter 
of direct verification to prove that the asymptotic behavior of the function S{E,h) is C2{h)/E'^ with C2{h) an E- 
independent function. Hence it follows that 5{Ey,h) — (1 — v)^5{E,h) oc (1 — v)'^vE~^ as i? — > 00. Thus, putting 
6^^\E, /i) = as a first approximation for the function 5{E, h), the second one can be found from the equation 



l-BiE)^-R2iE,k) 



5^'^\e, h) =^i{E, h), 



where we introduced 



Ri{E,h) = Y^j <^k{v,E,){n{E,E,,h){l~vy -[l+uj{E,h)v]]r^{v,E)dv, l>0 
k J° 

and ^i{E, h) = Ro{E, h). Repeating the consideration, one can proof by induction that S{E, h) — S^'''^ {E , h) cx ci{h)/E^ 
as E 00. Let us define 

ei{E,h) ==(5(')(£;,/i)-(5('-i)(£;,/i), i>2, 

^i{E,h)^y] [ <S>k{v,E,)n{E,E,,h)[ei{E,,h)-(l-vyQi{E,hMv,E)dv, I > 2. 
k Jo 

Then the following recursion chain of equations for the functions Qi{E, h) is derivable from the above reasoning: 

— -B{E)—~Ri{E,h) Qi{E,h)^%-i{E,h), l>2, (6.5) 
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The solution to Eq. ( |6.5| ) is given by 



ei{E,h) 



exp 



Ri {£o{E,h^ h"),h")dh" 



^i-i{£aiE,h- h'),h')dh', 



where £o{E, h) is the root of Eq. ( p!^ ) with = in its left-hand side. 

The formal convergence of this procedure can be proved under quite general assumptions on the energy dependence 
of the functions involved; specifically if the functions B{E) and $/c(w, E) increase monotonically and sufficiently slowly, 
while 'Dq(E) decreases with energy so that CC-E") is a slightly varying function of energy. It follows from our numerical 
analysis that in "real environment" the rate of convergence is very high: the first approximation {S{E, h) = 0) works 
with a reasonable exactness up to about 6 km w.e. and 3-4 iterations are suffice to obtain a few-percent accuracy for 
the differential muon spectrum at /i < 18 km w.e. and E > 1 GeV. 

The results obtained by the method being discussed agree well with our previous calculations |126]. At /i < 16 km 
w.e. of standard rock, the m etho d was verified by the direct Monte Carlo calculation, using an updated version of 
the code by Takahashi et al. |125]. The accuracy of our calculation for the muon depth-intensity relation (DIR), 



V^{E,h)dE (Ea^^lGeV), 



is estimated to be under (2-3) % at all depths of interest. The systematic errors can only be caused by the uncertainties 
in the input parameters, namely, the boundary muon spectrum and the muon-matter interaction cross sections. An 
additional error arising on the comparison with the data of a particular experiment, is related to the u ncert ainties in 
the averaged density and chemical composition of the matter overburden ((p), (Z), (A), (Z/A), (Z'^/A)) |129|. Strictly 
speaking, the approximation of homogeneous medium may also introduce a systematic error into the calculation for 
the real inhomogeneous media |l3C 



VII. CALCULATED MUON DIR VS UNDERGROUND AND UNDERWATER DATA 

A. Early underground experiments 



In Figure we present a comparison between the calculated vertical intensity (vs. depth underground) of con- 
ventional muons and the data obtained in early underground experiments performed with relatively small detec- 
tors pl|-|25 ,p8|-^ as well as the Crouch's 1987 "World Survey" data jSOj. To expand the comparison, we represent 
in Figure]? a fragment of the same information relevant to shallow depths. 

The data obtained by Wilson and by Clay and Van Gemert in the late 1930s are rather uncertain since 
the techniques used were unable to estimate the effects of showers, scattering and 5-electrons. We have normalized 
these data to our curve. All the other data points in Figs. ^ and are absolute. The Crouch World Survey comprises 
the data of different experiments, in particular, the early KGF data and extensive data from East Rand 

Proprietary Mine (ERPM) near Johannesburg |3J]) at great depths (all the points at ft. > 7.5 km w.e.). All the data 
were converted by Crouch to standard rock {Z = 11, A = 22, p = 2.65 g/cm'^) with some correction for the depths. 
Crouch's original compilation also includes the data from the depths well beyond 18 km w.e., where the atmospheric 
muon contribution is entirely negligible compared to the neutrino induced muon flux (see below), as well as three 
data points from an underwater experiment |5J] (we dropped these three points intending to discuss the complete set 
of underwater data below). 

At /i > 11 km w.e., the flux 2^ of muons produced by atmospheric neutrino interactions in the surrounding rock 
becomes important. The value of can significantly vary from one experiment to another due to different registration 
thresholds, the topology of the matter overburden, and so on. So, to account for the neutrino-induced background, 
we shall use the specific experimental data rather than some theoretical predictions. In Figure ^ we use the result of 
Ref. [|0): I^, = (2.17 ±0.21) x lO'^^ cm-^s-^sr-i. 

According to Crouch, the presented data at ft > 1 km w.e. can be approximated by the following empirical function: 



I^(ft) = exp(Ai + A2h) + exp(A3 + A^h) 



(7.1) 

-0.001213 ± 0.000021 (the result of 
a least square fit). Here Ip(ft) is in cm^^s^^sr^^ and ft (the depth in standard rock) is in hg/cm^ (1 hg/cm^ = 1 m 



with A\ 



11.22 ±0.17, A2 



-0.00262 ±0.00013, A3 



-14.10 ±0.14, Aa 



w.e.). The fit (7.1) is in good agreement with the result of the Utah group |32| 
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FIG. 7. The same as in Figure o but for shallow depths. 
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As Figs. |6|,|7| suggest, our DIR for conventional muons agrees well with most of the data within a wide depth range 
from about 30 to 9000 m w.e. (the exceptions are the points from Refs. and at /i « 300 m w.e. and also the 
points from Ref. lying in the range 2.1 to 2.5 km w.e.). The maximum disagreement with the best fit (7.1) at 
/i = (1 ^ 7.5) km w.e. is about 10 %. However, a,t h — (9.5 12) km w.e., our intensity noticeably exceeds the ERPM 
data. At h = 11.5 km w.e., the disagreement ranges up to about 77% (or about 30%, if one take the experimental 
errors into account). Such an error goes far beyond the expected accuracy of our calculations and thus are attributable 
either to uncertainties in the input parameters (primary spectrum?) or to some systematics in the ERPM data. It 
is clear that the data give no indication of some PM fraction in the muon flux so we do not show the corresponding 
curves here. As an illustration, let us note that, at /i = 12 km w.e., the calculated muon intensities with the PM 
contribution which result from the QGSM, RQPM, and VFGS are respectively 1.7, 2.3, and 3.3 times larger than the 
Crouch best fit. 



B. Kolar Gold Fields 

Figure ^ shows a comparison with the data obtained from several detectors located at different levels in the deep 
mine of the Kolar Gold Fields, Mysore State, South India (vertical telescopes at 745, 1500 and 3375 m w.e., a 
horizontal telescope at 3375 m w.e. and proton decay detectors at 6045 and 7000 m w.e.). In Ref. js^, the neutrino- 
induced background has been subtracted from the data using the measured angular distribution of muons. The four 
curves in Figure ^ represent our predictions for the muon intensities without and with adding the PM contribution 
from the QGSM, RQPM, and VFGS. Our calculations are done for the Kolar rock with {Z) = 12.9, {A) 26.9, 
{Z/A) = 0.495, (Z^/A) = 6.31, and (p) = 3.05 g/cm^. 

Up to 6-7 km w.e., one can see an excellent agreement between our predictions and the KGF data, irrespective of 
the PM flux model. Contrary to the data presented in Figure ^ the KGF muon DIR visibly exceed the calculated 
TT, if- muon intensity at ft. > 7 km w.e., hinting at some PM contribution. Both the RQPM and the VFGS model are 
in agreement with the KGF data up to about 10 km w.e^ but the VFGS model better fits the deeper data. This is not 
in contradiction with the situation presented in Figure |5| (b) for the sea- level integral spectrum when the ambiguities 
of the conversion procedure mentioned in Section [V| are taken into account. 



C. Baksan 

In Figure ^ we show a comparison of the calculation and the data obtained with the Baksan underground scintillation 
telescope (North Caucasus, Russia). The data obtained at zenith angles 50° - 70° (Ref. ^) and 70° - 85° (Ref. [|7)) 
were converted by the authors to vertical direction and to standard rock and the neutrino-induced muon background 
was subtracted from the data at high depths. A systematic difference between the two sets of data takes place in the 
depth interval from 6 to 9 km w.e.: in the first set (50° — 70°) a bump of intensity is clearly visible, while there is no 
such bump in the second set of data. 

The authors of Ref. [ p6| argue that the observed bump can be interpreted in terms of prompt muons. In our view 
this is not the case. The odds are that the bump is caused by errors in the determination of the oblique depths. 
Beyond the interval 6-9 km w.e., the data of both sets fall on a smooth curve. At the same time, the data at ft, > 10 
km w.e. may be attributed to the presence of some PM fraction in the measured underground muon intensity. Because 
of rather large experimental errors any model of charm production under consideration cannot be excluded by the 
Baksan data (including the case with the zero PM contribution), but it seems the data are more favorable for the 
RQPM. A collation of Figures |9| and || (b) suggests that the sea-level integral spectrum reconstructed in Ref. |^ from 
the Baksan muon DIR was distinctly overestimated. 



D. Mont Blanc Lab 

Figure |lO| shows the comparison of the predicted DIR for the conventional muons with the single muon intensity 
measured with the detectors SCE and NUSEX located in the Mont Blanc Laboratory. Our calculation rep- 

resents the muon intensity averaged upon the muon multiplicities and therefore we can make nothing more than 
qualitative conclusions from the comparison. 

In the overlapping region {h < 7 km w.e.) the data of both detectors superimpose and (with allowance for the 
multi-muon events) agree with our result. However, as is clear from the figure, the NUSEX DIR has a much greater 
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FIG. 8. Muon intensity vs Kolar rock thickness. The KGF data are from Ref. jS^]. The curves are for the ■R,K-m\ion DIR 
and for the DIR with the PM contributions calculated according to the RQPM and VFGS. 
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FIG. 9. Muon intensity vs standard rock thickness measured with the Baksan scintillation telescope [ p6|j3^ ]. The notation 
for the curves are the same as in Figure ^. 
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FIG. 10. Single muon intensity vs standard rock thickness measured in two experiments under Mont Blanc, SCE |Q and 
NUSEX ||l|. The curve is for the predicted Tr.if-muon DIR. 
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abrupt grade compared with the predicted one for the tt, ii'-muon DIR and, at (9-^11) km w.e., the predicted intensity 
(without any PM contribution) is 2-3 times higher than the NUSEX data. So large a discrepancy has no relation to 
the multi-muon events, whose intensity decreases with depth quicker than the single-muon one and is negligible at 
h > 6 — 7 km w.e. within a few percent accuracy. 

It is our opinion that the NUSEX result at large depths is incorrect. Notice that the muon DIR measured in the 
NUSEX experiment has been converted to standard rock. Although the averaged values of p, Z, A, Z/A and Z'^ /A in 
the Mont Blanc rock are rather close to the "standard" ones, this conversion might be a serious source of a systematic 
error because of very complicated and heterogeneous (layered) chemical composition of the rock (see Ref. [p5|). 

We note here that our calculations are in good agreement with the result of the French-American muon experi- 
ment ^] also carried out in the Mt. Blanc tunnel with a GM telescope. The depth range explored in that experiment 
was 0.5 to 5 km w.e. and therefore overlaps in part the SCE-NUSEX depth range. This suggests that the SCE and 
NUSEX experiments may have an added source of systematics related to their experimental procedures. At the same 
time, the recent NUSEX measurements of the averaged muon energy underground [131| are in good agreement with 
our predictions. 



E. SOUDAN 1/2 

Figure |ll| represents the comparison of our prediction with the data from SOUDAN 1 and SOUDAN 2 underground 
experiments [^,^ (the data points are taken from the compilation presented in Ref. |Q). The SOUDAN data were 
normalized to DIR for standard rock using the Crouch World Survey, as described in Ref. . 

Despite some spread of experimental points and a bump at 3 -f- 4 km w.e., one can see a reasonable agreement 
between the calculated tt, K-mnon DIR and the data up to about 7 km w.e., but the last point (^ 8.4 km w.e.) is 
almost 2.5 times below the predicted curve, as in the case of the NUSEX data. 



F. Frejus 

In Figure |l2|, we compare our calculations with the data of the Frejus detector |3^ (the underground laboratory 
was located in a tunnel of the same name under the Alps). The Alpine rock thickness has been converted into hg/cm^ 
of standard rock. We do not include in our assemblage the new and very detailed data from the Frejus detector 
recently reanalyzed in Ref. (see also Ref. [^). The point is that the original data sample has been subdivided 
into throughgoing, multiple and stopping muons. These subsamples are very dependent of the features peculiar to 
the experiment and thus cannot be directly compared with our calculations. 

According to Ref. pO[ |, the neutrino-induced muon background becomes dominant at /i > 13 km w.e. and the 
measured mean background flux is — (3.67 ± 0.66) x 10^^'^cm^^s~^sr~^. One can see that the intensity calculated 
without the PM contribution and corrected for this background fits the Frejus data almost everywhere, although in 
the vicinity of /i = 10 km w.e. some hint of an excess over the data (similar to the more evident excess in ERPM 
and NUSEX) does take place. The PM contributions calculated with RQPM and VFGS do not fall into Frejus data. 
However, in view of the experimental uncertainties there is no telling that the RQPM prediction is in serious conflict 
with the Frejus result. Clearly the same is all the more true for the QGSM. 



G. Gran Sasso Lab 

The recent data from the two largest underground detectors MACRO [Q and LVD (located in the Gran Sasso 
Laboratory) are presented in Figs. |l^ and ^ respectively. The data of MACRO are converted to standard rock. 
The error bars include statistical uncertainty, systematic uncertainty for the topographical map and the additional 
estimated systematic scale uncertainty of ±8 %. Taken alone, the statistical errors in the MACRO experiment are 
very small. The main contribution to the absolute scale uncertainty comes from the assumption of a homogeneous 
mountain instead of a layered structure. The data of LVD are presented both for the Gran Sasso rock and standard 
rock. Errors include both statistical and systematic uncertainties. Notice that at h < 5 km w.e., the statistical errors 
are less than the size of the circles in the figure. 

The depths currently accessible for observation with the detector MACRO are insufficient to study prompt muons. 
Thus, in Figure |l^ we present the calculated tt, isT-muon DIR alone. In contrast, the LVD data (Figure |l^) overlap 
the total depth range where the PM contribution might be essential. We use for our calculated curves the following 
value of the neutrino-induced muon background: = (2.98 ± 1.15) x 10~^^ cm~^s~^sr~^ p5[ |. 
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FIG. 11. Muon intensity vs standard rock thickness measured in the SOUDAN 1 and SOUDAN 2 experiments [^2[^. The 
curve is for the predicted tt, if-muon DIR. 
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FIG. 12. Muon intensity vs standard rock thickness measured in the Frejus experiment [^. The dashed curve represents our 
TT, _R'-muon DIR, the solid curve represents the same plus the neutrino- induced muon background, , according to Ref. . 
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FIG. 13. Muon intensity vs standard rock thickness measured in the experiment MACRO |^^. The curve is for the calculated 
TT, K-muon DIR. 
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FIG. 14. Muon intensity vs standard rock and Gran Sasso rock thickness The notation of the curves is the same as in 
Figure |l^ but with the neutrino- induced muon contribution deduced in Ref. Q . 



34 



As Figure 13 suggests, in the range from 3200 to about 6000 hg/cm^ the MACRO data are systematically in excess 
of our predicted muon intensity by about (8-10) %, what is beyond the total systematic error estimated in Ref. Q. 
This fact seems to be in dramatic contradiction with the sea-level muon spectrum reconstructed from the MACRO 
underground data under discussion (see Figure Indeed, our sea-level spectrum of the tt, muons is in good 
agreement with the MACRO fit (5.1) of the sea- level spectrum from 600-700 GeV up to 6-7 TeV and it stands ou t 
above the MACRO fit at higher energies. However, as noted in Section 0, the overall systematic error of the fit (5.1) 
is large enough to explain this contradiction, at least formally. 

At all depths, the LVD data are in excellent agreement with our calculations for the conventional muon DIR. 
Therefore, the data favor the models of charm production which predict a very low PM contribution (QGSM, pQCD, 
DPM). However, the RQPM cannot be excluded by the LVD result as yet. These conclusions concur with the 
conclusions of Ref. . The consistency of the data with our calculations for standard and Gran Sasso rocks provides 
an important confirmation for the correctness of the new conversion procedure to standard rock used in the LVD 
analysis BtI]. 



H. Underwater data 



Some problems of the underground muon experiments can be overcome by measurements underwater (and "un- 
derice") owing to unlimited (in principle) detection volume, uniformity and well known composition of the matter 
overburden. 

We present the total (to our knowledge) assemblage of underwater data in Fig. [l5[ The measurements with 
compact closed installations were performed in Suruga-bay, West Pacific (Higashi et al. [^l[), in Lake Geneva (Rogers 
and Tristam [^ ), in the Atlantic Ocean, Black, Mediterranean, and Caribbean Seas during several expeditions 
of research ships (Davitaev et al. j52j, Fyodorov et al. Q]). The measurements with open detectors (strings with 
phototubes), the prototypes of future large-scale neutrino telescopes, were performed in the Pacific Ocean off the West 
coast of the island of Hawaii in 1987 (the DUMAND Short Prototype String, Babson et al. [^), in the Mediterranean 
Sea a short way off Pylos, during three expeditions in 1989, 1991 and 1992 (the NESTOR prototypes, Anassontzis et 
al. [|56[|), in Lake Baikal during two expeditions in 1992 and 1993 (the stationary prototypes of the underwater neutrino 
telescope NT-200, Belolaptikov et al. ||5^,||), and at the South Pole with AMANDA (AMANDA-B4 experiment ||9|). 

Our calculation for the vr, if-muon DIR was done for sea water with (Z) = 7.468, (A) = 14.87, (Z/A) = 0.5525, 
(Z^/A) = 3.770 and (p) = 1.027 g/cm^. At h<7 km, the difference with the DIR for pure H2O is less than 1 % and 
can be neglected as compared to the theoretical and experimental uncertainties. There are two predictions in Fig. |l^: 
upper curve corresponds to muon threshold energy of 1 GeV and lower one corresponds to 20 GeV. 

At shallow depths (to 175 m) there are two measurements with very good statistics (Higashi et al. |^,^, Rogers 
and Tristam |Q), but the results of Higashi et al. are (except for the inclined data points at 105 m) lower by 15 to 
30 % than the result of Rogers and Tristam. According to Ref. ||5^ , one reason for the discrepancy is believed to be 
as follows. Higashi et al. normalized their data to an intensity derived from earlier underground measurements and 
measurements of the sea-level muon spectrum. The intensity chosen for the normalization is not quoted in Ref. , 
but was almost certainly too low. Our prediction is in excellent agreement with the absolute intensity obtained by 
Rogers and Tristam. This provides good support of our nuclear cascade model at low energies. However, the absolute 
measurements of Davitaev et al. |5^] are systematically lower than our prediction at < 1 km. 

As for greater depths, (1 4) km, it can be concluded that our prediction is in tolerable agreement with the data 
from the DUMAND and NESTOR prototypes as well as with the data of Fyodorov et al.; the discrepancy with a few 
specific data points is within (1 — 1.5)(t and is compatible with the overall data scattering. The data of the Baikal 
Collaboration [|7|j58) and the AMANDA Collaboration |5^ are in very good agreement with our curve. 

As is evident from the foregoing, the present-day state of the large-scale underwater projects does not permit 
to compete with the underground detectors as yet. In particular, the (slant) depths explored by the present-day 
underwater experiments are too small to get useful information on the PM flux. It is hoped that the situation will 
change in the immediate future. 
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VIII. CONCLUSIONS 



In this work we have attempted to study the vertical flux of high-energy cosmic-ray muons "from top to bottom" , 
that is from the primary spectrum of cosmic rays to underground/water muon intensity. Based upon the comparison 
of our calculations with the present-day ground-level, underground and underwater measurements, we have reached 
the following conclusions. 

Below 1 TeV/c, the spread of the data on the vertical sea-level muon spectra (differential and integral) measured 
in different experiments runs up to about 50 % and it is as much as (25-30) % even among the data of absolute 
measurements. Our calculations in this momentum range most closely fit the absolute 1984 data from the Nottingham 
magnetic spectrograph jlj] . Below 5-6 TeV/c, they agree with the indirect data deduced in most of the underground 
experiments (Artyomovsk Baksan |37|6J], MACRO jii), Frejus [Q, MSU |Q) without reference to a charm 
production model. 

All available indirect sea level data (Baksan, KGF, Frejus, MSU) show flattering of the sea-level muon spectrum 
at energies above ~ 20 TeV. The basic conclusion of Refs. [ ^ , ^8p^j6^ is that this flattering is due to the charm 
production in the atmosphere. This conclusion can be deduced from the experimental and theoretical results summa- 
rized in Figures |b and |b. First one should note that almost all experimental works use a power law muon spectrum 
with the spectral index similar (more or less) to our one at low energies. Further, they assume this index energy 
independent for conventional muons. According to these assumptions, a prompt muon contribution is the main cause 
of the flattering. 

However, the analysis shows definitely that the sea level data in the aggregate cannot be quantitatively described 
by a single charm production model. What this means is that the results of different experimental groups are in 
rather poor agreement with one another. Essentially all data on the differential spectrum shown in Figure |b lead to 
the conclusion about very high charm production rate, as is predicted by the VFGS model or even higher. At the 
same time, the data on the muon integral spectrum from different groups require different rates of charm production: 
the KGF data Q do not contradict to the QGSM or RQPM, while the Baksan data clearly favor the VFGS model. 

At the depths from about 30 m w.e. up to 6-7 km w.e., essentially all underground data on the muon DIR correlate 
with each other and with the predicted intensity for conventional (tt, if) muons, to within 10%. Hence it follows 
that our nuclear cascade model is valid with the same precision from about 8 GeV up to 4-5 TeV of muon energy 
at sea level, i.e. up to about 100 TeV/nucleon in the primary spectrum. This precision is distinctly better than 
one might expect with the current uncertainties in the input parameters (including the primary spectrum model and 
the muon-matter interaction cross sections). It is important that the underground data at ft, < 6 — 7 km w.e. do 
more than correlate well, but also have a very good statistical accuracy. Therefore, they may be of utility, among 
other things, for a normalization of the atmospheric neutrino flux in the range of neutrino energies most essential for 
j/-induced muons. 

The present-day world underwater results, though moderately detailed, provide a very important check upon the 
accuracy of the underground experiments, since they are free of the uncertainties in the density and composition of the 
matter overburden. The data obtained with the prototypes of future large underwater neutrino telescopes, especially 
with the Baikal Neutrino Telescope and with AMANDA, are in good agreement with the underground data and with 
the present calculations. 

The situations with the underground data at ft > 7 km w.e. is unsatisfactory in the same sense that it is with the 
ground-level data at high energies. The data from KGF and also from Baksan (measurements at i? = 50° — 70°) 
demonstrate clear excess over the predicted tt, if-muon curve providing fair indication of PM production. However, 
our calculations show that the KGF and Baksan sea-level data (Figures |^, |^) are inconsistent with their own (source) 
depth-intensity relations (see Figures ||, ||). Namely, the Baksan DIR is well explained by the RQPM rather than the 
VFGS model and the KGF depth-intensity curve is rather well described by the VFGS model. On the other hand, 
the corresponding sea-level spectra (Figure ^) show appreciably higher charm production rate in the case of Baksan 
and lower one in the case of KGF. Probably, these inconsistencies arise due to some uncertainties in the conversion 
procedures. 

The data points from ERPM and especially from NUSEX |Q are under the tt, if-muon curve. The data of 

all the other underground experiments are between the extremes represented by the Baksan/KGF and NUSEX. In 
particular, the most recent LVD data p7| are in good agreement with our tt, /C-muon DIR and therefore they favor 
a low charm production rate (as predicted by the QGSM or lower). 

In closing, no undisputed conclusion about PM production can be extracted from the world underground data. 
Considering that the statistical significances of the underground results at great depths are comparable and quite 
tolerable, this situation is attributable to the fact that certain of the experiments discussed have unrecorded systematic 
uncertainties. 
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APPENDIX A: SPECTRA OF MUONS FROM INCLUSIVE DECAY OF D AND Ac 

1. -D — » nVfiX decay 



1(1 - rl) (1 - 5rl - 2r%) + r^^x - i(l - 2rl)x^ + ^x' 



rjj In 



' D 



Zd = ^{1- r%) (1 - 8r|, + r%) - r% Inrf,, 

Here r = s'^/mjj and s^^' is the effective invariant mass square in the decay. The best fit to the data on the 
differential and total decay rates is achieved using 

0.63 GeV for D+ and £>" 
^ 0.67 GeV for D° and D° 

Therefore rn± ~ 0.337 and r^o^^o ~ 0.359. 

2. Ac — > fJiVfiX decay 



l<«<j<3 



l<i<J<3 
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sen{x) 



se22{x) 



3S23(a;) 
3333(2;) 



11 



21 



y - 3rA - yrl - 12ri 



9/ 



9r' 



-rl + 6ri(2 + rA)x - - (3 - 2rA - Sr^) x'^ 



(3 + 6rA + rj) 



+ri(l2 + 9rA + 6ri-ri)ln 



{1-x 



eei2{x) = 



aei3(a;) 



6rA - 9ri - 24ri + 3ri + 18ri - -rl + 2Arlx - 3 (l - 2rA - ri) x'' 



2ri (3 + 6rA + ri) 2ri 



1 -a; 



{1-xy 



+ 2ri(l2 + 3rA + 6ri-ri) In 



+ 3ri - 3rl - ^rl - Urix + 3(1 - 3rl)x^ - ^x^ 



+ 



2r-i(3 + ri) 



2rl 



2ri(3 + ri)ln 



1-x 



1 



(l + rA)2 



23 



. 40 fi 22 7 4 o 

24ri - 21ri + -rl + -rl - -ri 



+2rl (6 + 9rA + 6r 



rl)x~-{l- 4ri 



4rl+rt)x^ 



^(1+rA-ri); 



(3+12rA + 14ri + 8ri + r4) r6(l + rA)^ 



+rl (12 + 21rA + 24ri + lOri + Ari - rl) In 
11 



(1-x, 
1 - a; 



^ + 3rA - yri + 12ri + ^ri - 9rl + ^-rl - 6ri(2 - rA)x - ^ (3 + 2rA - 5ri) x^ 



I 83 1 ^A(3-6rA + ri) 
^3" 1-ar (1-^)2 



- ri (12 - 9rA + 6ri 



In 



1 



511 = (1 - rl) (1 - 2rA - 7rl - 20ri - 7rl - 2rl + rl) ~ 2Arl (l + ri + rl) InrA, 

512 = (1 - ri) (1 - 4rA - 7ri - 40ri - 7ri - 4ri + rl) - 24rl (2 + ta + 2ri) InrA, 
aei3 = ■^23 = 0, 



3322 



1 



(1 + rA)2 



- - rA - 6ri - 28ri - 32ri + 28rX + 6ri + 
5 



--r 



10 



24ri (1 + 2rA + 3ri + 2ri + rl) InrA 



3833 = (1 - ri) (1 + 2rA - 7ri + 20ri - 7ri + 2ri + rl) + 24ri (l - rA - ri) InrA- 

Here ri = s^f/'^Ac ^^'^ ■^x ^he effective invariant mass square. The best fit to the data on the differential and 
total decay rates is achieved using ^J's^ = 1-27 GeV. Therefore rA ~ 0.551. For the decay form factors averaged over 
^2 we have /i w 0.991, /2 w 2.170, and /s w 0.805. 
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APPENDIX B: MUON-MATTER INTERACTIONS AT HIGH ENERGIES 



Here we present with some comments a listing of the cross sections for the muon-matter interactions and the 
formula for ionization loss. In what follows, Z and A are the atomic number and atomic weight of the target nucleus; 
E and E' are the respective energies of initial and final muons; v is the fraction of energy lost [E' = E{\ — v)); a and 
re are the fine structure constant and the classical electron radius; nie and are the electron and muon masses; 
c= 1. 



1. Direct e^e pair production 



The direct p air p roduction cross section goes roughly as to over most of the range, v > 0.002 (see, for 
example, Ref. [ 132 |). Because of this, the energy loss through pair production is usually considered as continuous. 
Nevertheless, as it follows from our calculations, the fluctuation effect related to this process is not negligible and 
it grows in magnitude with depth. In the present work, only a part of the direct pair production cross section, 
corresponding to relatively large energy losses (v > vq = 2 x 10""*), is included into the collision integral of the muon 
transport equation while the energy losses caused by the range v < vq were treated as continuous. The exact (in the 
ultrarelativistic limit) numerical results for dup/dv have been obtained by Kel'ner and Kotov and can be found in 
Ref. WOl. In this paper, we use the following simple approximation of the exact results: 



^do, ^ 16 2 
dv TT 



^ , 1.7 X l{)-^(v + 1.05 X 10-*) 
v(v + 0.006)^ 



cxp [-0.025 ln2(£:/m^)] 



0.323 ln(103u) 



•f{E,v), 



where 



f{E,v) 



l + 2v 
l + k/{vE)'' 



k = 0.02 GeV, 



for 10^"^ <v< 0.2 and f{E,v) = 1, for vq < v < 10^^ or v > 0.2. A weak (logarithmic) Z-dependence of the exact 
function F{E,v) is neglected in this approximation. The accuracy of the approximation is better than 10% within 
the most essential interval of v {vq < v < 0.1). It provides a reasonable (a few-percent) precision for the calculated 
muon DIR. 



2. Bremsstrahlung 



We use the formula derived by Andreev et al. [133| with regard to the nuclear target structure and the exact 
contribution to the cross section given by atomic electrons: 



dab 
dv 



~ a. 2r„Z- 



{2-2v + V^) *i (g^i„, Z)--{l-v) ^2 (gmin, Z) 



*1,2 (gmin, Z) = -^^2 ilmin, Z) - Ai^2 ((Zmin, Z) . 



*0(g„,„,Z) 1 + ln^ 



*2 (ftnin, Z) 




1 1 

— xi arctan 1 

1 / xi Z 



1 / , fnla-l \ 1 

— 1 + In — - — - — X2 arctan — , 

2 \ 1 + X2 I X2 



2xi { 1 — xi arctan ■ 



3 xi 
xi 4 1 + xf 

2xi I 1 — X2 arctan h - In — 

a;2 4 1 - 



.^2 
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A ^ 7\ 1 "^A* , C, C+1 

Ai (q,nin, Z) = In h - In -, 

Qc 2 C - 1 



A2 (qmin-.Z) = In- 



3-C In- 



mf.v 



2£;(l-w 



a-i 



111.7 



a2 



724.2 



1.9to^ 



From the foregoing equations it can be shown that, in the hmit of complete screening, i.e. for 



lz{v,E) = 



200(?n 



1/3 



ITeV 



E / l-t; 



< 1 



(where 7^ is the degree of screening), the bremsstrahlung cross section is a function of the variable v only (scaling). 
However for values of v which are not too small (namely, at 1 — w <C 1) complete screening occurs only at very high 
energies, E ^ IQ TeV. At lower energies, the cross section grows logarithmically with E. It should be noted that the 
same estimate [E ^ 10 TeV) is also true as a limit of full screening for the pair production cross secti on. 

Now, let us discuss briefly the corrections to the Born approximation. Recently, it was shown |134] that there are 
two such corrections: in the region of small (g ~ rn^/ E) and large {q ~ m^^) momentum transfers which correspond 
to large and small impact parameters, respectively. The correction from the first region (large q) is just t he s ame as 
in the case of electron bremss trah lung (it is the well-known correction of Davies, Bethe, and Maximon [135|). The 
essentially new result of Ref. |134] is that the second correction has the opposite sign and nearly compensates the 
first one. As a result, the Born approximation formulas have rather good accuracy for the muon bremsstrahlung even 
for very heavy targets. In the case of interest for the present study, the corrections under discussion prove to be 
completely negligible. 



3. Photonuclear Interaction 



For the photonuclear interaction of muon we use the generalized vector dominance model (GVDM) |136|. Within 
the vector meson dominance hypothesis, the differential cross section for muon photonuclear interaction, dan/dv, 
is expressed in terms of the total cross section for virtual photon absorption by nucleons and nuclei. The GVDM 
adequately describes the features of these cross sections in the diffraction region (low 4-momentum transfers, Q^, 
and large photon energies, v): a growth with energy o f the cross section for nucleon photoabsorption and shadowing 
effects in nuclear photoabsorption. According to Ref. [136|, 



dv 



Stt 



Aajp{iy)v < H{v) In 1 + 



2ml 



In 













mlJ_ 



G{z) 



H{v) In 1 + - 



t J m\+t 



2ml 



1 - 



0.25ml 



9 



Here v = vE is the virtual photon energy and 



H{v) = 1 - - + 



G(z) = 



(1 + z)e-^ - 1 



= 0.00282A^/^c7-yp{iy 



t = 



9 2 



m\ = 0.54 GeV^ 



= 1.80 GeV^ 
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The differential cross section is proportional to the total cross section for absorption of a real photon of energy 
V — s/2m N = vE by a nucleon, <Tjn. In the present calculations, we adopt the Regge-type parametrization for cr^jv 
from Ref. p^ , 



c-yAT — [67.7s 



0.0808 



129s 



-0.452.51 



^b 



(Bl) 



(s in GeV^). This model give s the best fit to the accelerator data. At y/s - 
parametrization of Ref. [1_36| used in our previous calculations, 



200 GeV it differs by about 9 % from the 



cr.^jv - [114.3+ 1.647 ln^(0.0213z^)] 



(B2) 



in GeV). The disparity in the DIR resulting from the difference in the models (Bl) and (B2) for a-^N, is completely 
negligible up to ^ 10 km w.e. (independent of rock composition) and it is small at greater depths. Namely, the muon 
intensities calculated with the use of the parametrization (Bl) exceed those calculated with the parametrization ( p32| ) 
by 1.2, 2, 3, and 5% at respectively 12, 14, 16, and 18 km w.e. of standard rock what is of no importance for the 
interpretation of the current underground data. 

The growth of a-yN with the photon energy causes da„/dv to depend on the muon energy, E = v/v. The shadowing 
effect of nucleons inside a target nucleus gradually compensating the energy dependence of cr^jv, but a logarithmic 
growth of dcTn/ dv quantitatively remains up to £^ ~ 10 TeV and possibly in the asymptotics. 

One should note here that, within the VDM approach, the growth of u-^m with the photon energy is resulted by 
the growth of the hadron-nucleon cross section. However, non-VDM corrections to a-ff^ may, in principle, be not 
neghgible. A part of these corrections is caused by the pQCD ("minijet") contribution to the 77V total cross section 
being determined by the proton perturbative structure function (rather than an intermediate vector meson exchange) . 
This correction is small because the pQCD cross section is dominated by the "VDM photon" |138]. The second non- 
VDM correction is dominated by direct photon-proton reaction; it corresponds to the so-called "unresolved photon" . 
The magnitude of this correction strongly depends on the poorly known behavior of the gluon structure function in 
the proton, gp{x), at small x and, of course, on the usual QCD parameter p^'". For example, if gp{x) cx x""^'^, the 
correction from the 73 qq subprocess behaves with the photon energy as ^/s ]139|] and depends on p™'" as (p™'")_^. 

Available cosmic-ray data obtained with underground detectors |140| | (for ^ 10 TeV) and with EAS arrays |141] 
(up to 10'^— 10* TeV) is in agreement with formulas (Bl) or (B2). Nonetheless, the photonuclear interaction rests one 



of the sources of uncertainties at very high muon energies and, at the same time, a noteworthy subject for investigation 
with the future large-scale underwater telescopes. Luckily, this uncertainty is of little importance for the muon DIR. 



4. Ionization energy loss 



The ionization loss of a muon of energy E is given by the Bethe-Bloch stopping-power formula corrected to the 
density effect pi (see also Ref. pi), 



dE\ 
dx); 



CoZ 
/32 A 



In 



2meP^Wa 



AE^ 



20^ -5 -U 



Here Cq = 0.1535 MeVg ^ cm^, p is the muon momentum, l3 = p/E is the muon velocity. 



2rneP^ 



TO^ -I- TTlg -I- 2meE 



is the maximum energy transfer from the muon to an atomic electron. The function S is the density-effect correction. 
Its numerical values are given by the Sternheimcr's fit formula |142|, 



d = e{X - Xq) [4.6052X + a0{Xi - X){Xi - X)"' + C] , 

where is the step function {0{x) = at a; < and 6{x) = 1 at a: > 0), X = log (p/m ^). The values Xq, Xi, 
a and m depend on the substance (for the specific values we used the data of Ref. |143| with some modifications 
in the cases of Baksan and Kolar rocks); C ~ — [2lii(Iz /hvp) + 1], where Iz is the mean excitation energy, hvp = 
28.816-\/ pZ/A is the plasma energy (in eV) and p is the density of the medium (in g/cm'^). The shell-correction 
term, U — 2Ck/Z + 2Cl/Z -I- . . . , is generally neghgible for the energies at which the density-effect correction 5 is 
significant. 
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